LA Storage: diagonal matrix storage

src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs

@@ -162,6 +162,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
+        internal DenseColumnMajorMatrixStorage<Complex> Storage
+        {
+            get { return _storage; }
+        }
+
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>
@@ -217,7 +222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             var denseTarget = target as DenseMatrix;
             if (denseTarget != null)
             {
-                _storage.CopyTo(denseTarget._storage);
+                _storage.CopyTo(denseTarget.Storage);
                 return;
             }
 
@@ -361,7 +366,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </summary>
         public override void Clear()
         {
-            Array.Clear(_data, 0, _data.Length);
+            _storage.Clear();
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs

@@ -31,6 +31,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
     using System.Numerics;
     using Generic;
     using Properties;
+    using Storage;
     using Threading;
 
     /// <summary>
@@ -45,6 +46,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
     [Serializable]
     public class DiagonalMatrix : Matrix
     {
+        readonly SparseDiagonalMatrixStorage<Complex> _storage;
+
+        /// <summary>
+        /// Gets the matrix's data.
+        /// </summary>
+        /// <value>The matrix's data.</value>
+        readonly Complex[] _data;
+
          /// <summary>
         /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class. This matrix is square with a given size.
         /// </summary>
@@ -54,7 +63,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </exception>
         public DiagonalMatrix(int order) : base(order)
         {
-            Data = new Complex[order];
+            _storage = new SparseDiagonalMatrixStorage<Complex>(order, order);
+            _data = _storage.Data;
         }
 
         /// <summary>
@@ -68,11 +78,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </param>
         public DiagonalMatrix(int rows, int columns) : base(rows, columns)
         {
-            Data = new Complex[Math.Min(rows, columns)];
+            _storage = new SparseDiagonalMatrixStorage<Complex>(rows, columns);
+            _data = _storage.Data;
         }
 
         /// <summary>
-        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class with all entries set to a particular value.
+        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class with all diagonal entries set to a particular value.
         /// </summary>
         /// <param name="rows">
         /// The number of rows.
@@ -80,13 +91,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <param name="columns">
         /// The number of columns.
         /// </param>
-        /// <param name="value">The value which we assign to each element of the matrix.</param>
+        /// <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)
         {
-            Data = new Complex[Math.Min(rows, columns)];
-            for (var i = 0; i < Data.Length; i++)
+            _storage = new SparseDiagonalMatrixStorage<Complex>(rows, columns);
+            _data = _storage.Data;
+
+            for (var i = 0; i < _data.Length; i++)
             {
-                Data[i] = value;
+                _data[i] = value;
             }
         }
 
@@ -99,7 +112,8 @@ 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)
         {
-            Data = diagonalArray;
+            _storage = new SparseDiagonalMatrixStorage<Complex>(rows, columns, diagonalArray);
+            _data = _storage.Data;
         }
 
         /// <summary>
@@ -111,16 +125,16 @@ 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))
         {
-            var rows = array.GetLength(0);
-            var columns = array.GetLength(1);
+            _storage = new SparseDiagonalMatrixStorage<Complex>(array.GetLength(0), array.GetLength(1));
+            _data = _storage.Data;
 
-            for (var i = 0; i < rows; i++)
+            for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j < columns; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     if (i == j)
                     {
-                        Data[i] = array[i, j];
+                        _data[i] = array[i, j];
                     }
                     else if (((array[i, j].Real != 0.0) && !double.IsNaN(array[i, j].Real)) || ((array[i, j].Imaginary != 0.0) && !double.IsNaN(array[i, j].Imaginary)))
                     {
@@ -130,14 +144,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
         }
 
-        /// <summary>
-        /// Gets the matrix's data.
-        /// </summary>
-        /// <value>The matrix's data.</value>
-        internal Complex[] Data
+        internal SparseDiagonalMatrixStorage<Complex> Storage
         {
-            get;
-            private set;
+            get { return _storage; }
         }
 
         /// <summary>
@@ -156,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 row == column ? _data[row] : 0.0;
         }
 
         /// <summary>
@@ -178,7 +187,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         {
             if (row == column)
             {
-                Data[row] = value;
+                _data[row] = value;
             }
             else if (((value.Real != 0.0) && !double.IsNaN(value.Real)) || ((value.Imaginary != 0.0) && !double.IsNaN(value.Imaginary)))
             {
@@ -220,7 +229,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </summary>
         public override void Clear()
         {
-            Array.Clear(Data, 0, Data.Length);
+            _storage.Clear();
         }
 
         /// <summary>
@@ -247,13 +256,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 return true;
             }
 
-            if (diagonalMatrix.Data.Length != Data.Length)
+            if (diagonalMatrix._data.Length != _data.Length)
             {
                 return false;
             }
 
             // If all else fails, perform element wise comparison.
-            return !Data.Where((t, i) => t != diagonalMatrix.Data[i]).Any();
+            return !_data.Where((t, i) => t != diagonalMatrix._data[i]).Any();
         }
 
         /// <summary>
@@ -264,14 +273,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// </returns>
         public override int GetHashCode()
         {
-            var hashNum = Math.Min(Data.Length, 25);
+            var hashNum = Math.Min(_data.Length, 25);
             long hash = 0;
             for (var i = 0; i < hashNum; i++)
             {
 #if PORTABLE
-                hash ^= Precision.DoubleToInt64Bits(Data[i].GetHashCode());
+                hash ^= Precision.DoubleToInt64Bits(_data[i].GetHashCode());
 #else
-                hash ^= BitConverter.DoubleToInt64Bits(Data[i].GetHashCode());
+                hash ^= BitConverter.DoubleToInt64Bits(_data[i].GetHashCode());
 #endif
             }
 
@@ -351,7 +360,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
             else
             {
-                Control.LinearAlgebraProvider.AddArrays(Data, diagOther.Data, diagResult.Data);    
+                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);    
             }
         }
 
@@ -426,7 +435,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
             else
             {
-                Control.LinearAlgebraProvider.SubtractArrays(Data, diagOther.Data, diagResult.Data);
+                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
             }
         }
 
@@ -447,12 +456,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 throw new ArgumentNullException("source");
             }
 
-            if (source.Length != Data.Length)
+            if (source.Length != _data.Length)
             {
                 throw new ArgumentException(Resources.ArgumentArraysSameLength, "source");
             }
         
-            Array.Copy(source, Data, source.Length);
+            Array.Copy(source, _data, source.Length);
         }
 
         /// <summary>
@@ -474,12 +483,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 return;
             }
 
-            if (Data.Length != denseSource.Data.Length)
+            if (_data.Length != denseSource.Data.Length)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "source");
             }
 
-            Array.Copy(denseSource.Data, Data, denseSource.Data.Length);
+            Array.Copy(denseSource.Data, _data, denseSource.Data.Length);
         }
 
         /// <summary>
@@ -515,7 +524,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                     CopyTo(diagResult);
                 }
 
-                Control.LinearAlgebraProvider.ScaleArray(scalar, Data, diagResult.Data);
+                Control.LinearAlgebraProvider.ScaleArray(scalar, _data, diagResult._data);
             }
         }
 
@@ -559,12 +568,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
             else
             {
-                var thisDataCopy = new Complex[r.Data.Length];
-                var otherDataCopy = new Complex[r.Data.Length];
-                Array.Copy(Data, thisDataCopy, (r.Data.Length > Data.Length) ? Data.Length : r.Data.Length);
-                Array.Copy(m.Data, otherDataCopy, (r.Data.Length > m.Data.Length) ? m.Data.Length : r.Data.Length);
+                var thisDataCopy = new Complex[r._data.Length];
+                var otherDataCopy = new Complex[r._data.Length];
+                Array.Copy(_data, thisDataCopy, (r._data.Length > _data.Length) ? _data.Length : r._data.Length);
+                Array.Copy(m._data, otherDataCopy, (r._data.Length > m._data.Length) ? m._data.Length : r._data.Length);
 
-                Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, r.Data);
+                Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, r._data);
             }
         }
 
@@ -641,9 +650,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 result.Clear();
 
                 // Multiply the elements in the vector with the corresponding diagonal element in this.
-                for (var r = 0; r < Data.Length; r++)
+                for (var r = 0; r < _data.Length; r++)
                 {
-                    result[r] = Data[r] * rightSide[r];
+                    result[r] = _data[r] * rightSide[r];
                 }
             }
         }
@@ -691,9 +700,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 result.Clear();
 
                 // Multiply the elements in the vector with the corresponding diagonal element in this.
-                for (var r = 0; r < Data.Length; r++)
+                for (var r = 0; r < _data.Length; r++)
                 {
-                    result[r] = Data[r] * leftSide[r];
+                    result[r] = _data[r] * leftSide[r];
                 }
             }
         }
@@ -709,7 +718,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 throw new ArgumentException(Resources.ArgumentMatrixSquare);
             }
 
-            return Data.Aggregate(Complex.One, (current, t) => current * t);
+            return _data.Aggregate(Complex.One, (current, t) => current * t);
         }
 
         /// <summary>
@@ -722,7 +731,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         {
             // TODO: Should we return reference to array? In current implementation we return copy of array, so changes in DenseVector will
             // not influence onto diagonal elements
-            return new DenseVector((Complex[])Data.Clone());
+            return new DenseVector((Complex[])_data.Clone());
         }
 
         /// <summary>
@@ -790,24 +799,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         public override void CopyTo(Matrix<Complex> target)
         {
             var diagonalTarget = target as DiagonalMatrix;
-
-            if (diagonalTarget == null)
+            if (diagonalTarget != null)
             {
-                base.CopyTo(target);
+                _storage.CopyTo(diagonalTarget.Storage);
                 return;
             }
 
-            if (ReferenceEquals(this, target))
+            var denseTarget = target as DenseMatrix;
+            if (denseTarget != null)
             {
+                _storage.CopyTo(denseTarget.Storage);
                 return;
             }
 
-            if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, target);
-            }
-
-            Array.Copy(Data, diagonalTarget.Data, Data.Length);
+            base.CopyTo(target);
         }
 
         /// <summary>
@@ -817,7 +822,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         public override Matrix<Complex> Transpose()
         {
             var ret = new DiagonalMatrix(ColumnCount, RowCount);
-            Array.Copy(Data, ret.Data, Data.Length);
+            Array.Copy(_data, ret._data, _data.Length);
             return ret;
         }
 
@@ -871,9 +876,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
             // Clear the result and copy the diagonal entry.
             result.Clear();
-            if (columnIndex >= rowIndex && columnIndex < rowIndex + length && columnIndex < Data.Length)
+            if (columnIndex >= rowIndex && columnIndex < rowIndex + length && columnIndex < _data.Length)
             {
-                result[columnIndex - rowIndex] = Data[columnIndex];
+                result[columnIndex - rowIndex] = _data[columnIndex];
             }
         }
 
@@ -927,9 +932,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
 
             // Clear the result and copy the diagonal entry.
             result.Clear();
-            if (rowIndex >= columnIndex && rowIndex < columnIndex + length && rowIndex < Data.Length)
+            if (rowIndex >= columnIndex && rowIndex < columnIndex + length && rowIndex < _data.Length)
             {
-                result[rowIndex - columnIndex] = Data[rowIndex];
+                result[rowIndex - columnIndex] = _data[rowIndex];
             }
         }
 
@@ -937,21 +942,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// <returns>The L1 norm of the matrix.</returns>
         public override Complex L1Norm()
         {
-            return Data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
+            return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
         }
 
         /// <summary>Calculates the L2 norm.</summary>
         /// <returns>The L2 norm of the matrix.</returns>   
         public override Complex L2Norm()
         {
-            return Data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
+            return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
         }
 
         /// <summary>Calculates the Frobenius norm of this matrix.</summary>
         /// <returns>The Frobenius norm of this matrix.</returns>
         public override Complex FrobeniusNorm()
         {
-            var norm = Data.Sum(t => t.Magnitude * t.Magnitude);
+            var norm = _data.Sum(t => t.Magnitude * t.Magnitude);
             return Math.Sqrt(norm);
         }
 
@@ -968,7 +973,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         {
             var maxSv = double.NegativeInfinity;
             var minSv = double.PositiveInfinity;
-            foreach (var t in Data)
+            foreach (var t in _data)
             {
                 maxSv = Math.Max(maxSv, t.Magnitude);
                 minSv = Math.Min(minSv, t.Magnitude);
@@ -989,11 +994,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
 
             var inverse = (DiagonalMatrix)Clone();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                if (Data[i] != 0.0)
+                if (_data[i] != 0.0)
                 {
-                    inverse.Data[i] = 1.0 / Data[i];
+                    inverse._data[i] = 1.0 / _data[i];
                 }
                 else
                 {
@@ -1037,9 +1042,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
 
             result.Clear();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
         }
 
@@ -1102,9 +1107,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             }
 
             result.Clear();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
         }
 
@@ -1198,7 +1203,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 int end = Math.Min(columnCount, rowCount + columnInit);
                 for (var i = 0; columnInit + i < end; i++)
                 {
-                    result[i, columnInit + i] = Data[rowIndex + i];
+                    result[i, columnInit + i] = _data[rowIndex + i];
                 }
             }
             else if (rowIndex < columnIndex && rowIndex + rowCount > columnIndex)
@@ -1207,14 +1212,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
                 int end = Math.Min(columnCount + rowInit, rowCount);
                 for (var i = 0; rowInit + i < end; i++)
                 {
-                    result[rowInit + i, i] = Data[columnIndex + i];
+                    result[rowInit + i, i] = _data[columnIndex + i];
                 }
             }
             else
             {
                 for (var i = 0; i < Math.Min(columnCount, rowCount); i++)
                 {
-                    result[i, i] = Data[rowIndex + i];
+                    result[i, i] = _data[rowIndex + i];
                 }
             }
 
@@ -1228,9 +1233,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         public override Complex[,] ToArray()
         {
             var result = new Complex[RowCount, ColumnCount];
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             return result;
@@ -1379,9 +1384,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1447,9 +1452,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1511,9 +1516,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1568,7 +1573,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             var m = new DiagonalMatrix(order);
             for (var i = 0; i < order; i++)
             {
-                m.Data[i] = 1.0;
+                m._data[i] = 1.0;
             }
 
             return m;

src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs

@@ -162,6 +162,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
+        internal DenseColumnMajorMatrixStorage<Complex32> Storage
+        {
+            get { return _storage; }
+        }
+
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>
@@ -217,7 +222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             var denseTarget = target as DenseMatrix;
             if (denseTarget != null)
             {
-                _storage.CopyTo(denseTarget._storage);
+                _storage.CopyTo(denseTarget.Storage);
                 return;
             }
 
@@ -361,7 +366,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </summary>
         public override void Clear()
         {
-            Array.Clear(_data, 0, _data.Length);
+            _storage.Clear();
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs

@@ -31,6 +31,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
     using Generic;
     using Numerics;
     using Properties;
+    using Storage;
     using Threading;
 
     /// <summary>
@@ -45,6 +46,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
     [Serializable]
     public class DiagonalMatrix : Matrix
     {
+        readonly SparseDiagonalMatrixStorage<Complex32> _storage;
+
+        /// <summary>
+        /// Gets the matrix's data.
+        /// </summary>
+        /// <value>The matrix's data.</value>
+        readonly Complex32[] _data;
+
         /// <summary>
         /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class. This matrix is square with a given size.
         /// </summary>
@@ -55,7 +64,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public DiagonalMatrix(int order)
             : base(order)
         {
-            Data = new Complex32[order];
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(order, order);
+            _data = _storage.Data;
         }
 
         /// <summary>
@@ -70,11 +80,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public DiagonalMatrix(int rows, int columns)
             : base(rows, columns)
         {
-            Data = new Complex32[Math.Min(rows, columns)];
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(rows, columns);
+            _data = _storage.Data;
         }
 
         /// <summary>
-        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class with all entries set to a particular value.
+        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class with all diagonal entries set to a particular value.
         /// </summary>
         /// <param name="rows">
         /// The number of rows.
@@ -82,14 +93,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <param name="columns">
         /// The number of columns.
         /// </param>
-        /// <param name="value">The value which we assign to each element of the matrix.</param>
+        /// <param name="value">The value which we assign to each diagonal element of the matrix.</param>
         public DiagonalMatrix(int rows, int columns, Complex32 value)
             : base(rows, columns)
         {
-            Data = new Complex32[Math.Min(rows, columns)];
-            for (var i = 0; i < Data.Length; i++)
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(rows, columns);
+            _data = _storage.Data;
+
+            for (var i = 0; i < _data.Length; i++)
             {
-                Data[i] = value;
+                _data[i] = value;
             }
         }
 
@@ -103,7 +116,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public DiagonalMatrix(int rows, int columns, Complex32[] diagonalArray)
             : base(rows, columns)
         {
-            Data = diagonalArray;
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(rows, columns, diagonalArray);
+            _data = _storage.Data;
         }
 
         /// <summary>
@@ -116,16 +130,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public DiagonalMatrix(Complex32[,] array)
             : this(array.GetLength(0), array.GetLength(1))
         {
-            var rows = array.GetLength(0);
-            var columns = array.GetLength(1);
+            _storage = new SparseDiagonalMatrixStorage<Complex32>(array.GetLength(0), array.GetLength(1));
+            _data = _storage.Data;
 
-            for (var i = 0; i < rows; i++)
+            for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j < columns; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     if (i == j)
                     {
-                        Data[i] = array[i, j];
+                        _data[i] = array[i, j];
                     }
                     else if (((array[i, j].Real != 0.0) && !double.IsNaN(array[i, j].Real)) || ((array[i, j].Imaginary != 0.0) && !double.IsNaN(array[i, j].Imaginary)))
                     {
@@ -135,14 +149,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
         }
 
-        /// <summary>
-        /// Gets the matrix's data.
-        /// </summary>
-        /// <value>The matrix's data.</value>
-        internal Complex32[] Data
+        internal SparseDiagonalMatrixStorage<Complex32> Storage
         {
-            get;
-            private set;
+            get { return _storage; }
         }
 
         /// <summary>
@@ -161,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 row == column ? _data[row] : 0.0f;
         }
 
         /// <summary>
@@ -183,7 +192,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         {
             if (row == column)
             {
-                Data[row] = value;
+                _data[row] = value;
             }
             else if (((value.Real != 0.0) && !double.IsNaN(value.Real)) || ((value.Imaginary != 0.0) && !double.IsNaN(value.Imaginary)))
             {
@@ -225,7 +234,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </summary>
         public override void Clear()
         {
-            Array.Clear(Data, 0, Data.Length);
+            _storage.Clear();
         }
 
         /// <summary>
@@ -252,13 +261,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 return true;
             }
 
-            if (diagonalMatrix.Data.Length != Data.Length)
+            if (diagonalMatrix._data.Length != _data.Length)
             {
                 return false;
             }
 
             // If all else fails, perform element wise comparison.
-            return !Data.Where((t, i) => t != diagonalMatrix.Data[i]).Any();
+            return !_data.Where((t, i) => t != diagonalMatrix._data[i]).Any();
         }
 
         /// <summary>
@@ -269,14 +278,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// </returns>
         public override int GetHashCode()
         {
-            var hashNum = Math.Min(Data.Length, 25);
+            var hashNum = Math.Min(_data.Length, 25);
             long hash = 0;
             for (var i = 0; i < hashNum; i++)
             {
 #if PORTABLE
-                hash ^= Precision.DoubleToInt64Bits(Data[i].GetHashCode());
+                hash ^= Precision.DoubleToInt64Bits(_data[i].GetHashCode());
 #else
-                hash ^= BitConverter.DoubleToInt64Bits(Data[i].GetHashCode());
+                hash ^= BitConverter.DoubleToInt64Bits(_data[i].GetHashCode());
 #endif
             }
 
@@ -356,7 +365,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
             else
             {
-                Control.LinearAlgebraProvider.AddArrays(Data, diagOther.Data, diagResult.Data);
+                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
             }
         }
 
@@ -431,7 +440,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
             else
             {
-                Control.LinearAlgebraProvider.SubtractArrays(Data, diagOther.Data, diagResult.Data);
+                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
             }
         }
 
@@ -452,12 +461,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 throw new ArgumentNullException("source");
             }
 
-            if (source.Length != Data.Length)
+            if (source.Length != _data.Length)
             {
                 throw new ArgumentException(Resources.ArgumentArraysSameLength, "source");
             }
 
-            Array.Copy(source, Data, source.Length);
+            Array.Copy(source, _data, source.Length);
         }
 
         /// <summary>
@@ -479,12 +488,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 return;
             }
 
-            if (Data.Length != denseSource.Data.Length)
+            if (_data.Length != denseSource.Data.Length)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "source");
             }
 
-            Array.Copy(denseSource.Data, Data, denseSource.Data.Length);
+            Array.Copy(denseSource.Data, _data, denseSource.Data.Length);
         }
 
         /// <summary>
@@ -520,7 +529,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                     CopyTo(diagResult);
                 }
 
-                Control.LinearAlgebraProvider.ScaleArray(scalar, Data, diagResult.Data);
+                Control.LinearAlgebraProvider.ScaleArray(scalar, _data, diagResult._data);
             }
         }
 
@@ -564,12 +573,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
             else
             {
-                var thisDataCopy = new Complex32[r.Data.Length];
-                var otherDataCopy = new Complex32[r.Data.Length];
-                Array.Copy(Data, thisDataCopy, (r.Data.Length > Data.Length) ? Data.Length : r.Data.Length);
-                Array.Copy(m.Data, otherDataCopy, (r.Data.Length > m.Data.Length) ? m.Data.Length : r.Data.Length);
+                var thisDataCopy = new Complex32[r._data.Length];
+                var otherDataCopy = new Complex32[r._data.Length];
+                Array.Copy(_data, thisDataCopy, (r._data.Length > _data.Length) ? _data.Length : r._data.Length);
+                Array.Copy(m._data, otherDataCopy, (r._data.Length > m._data.Length) ? m._data.Length : r._data.Length);
 
-                Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, r.Data);
+                Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, r._data);
             }
         }
 
@@ -646,9 +655,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 result.Clear();
 
                 // Multiply the elements in the vector with the corresponding diagonal element in this.
-                for (var r = 0; r < Data.Length; r++)
+                for (var r = 0; r < _data.Length; r++)
                 {
-                    result[r] = Data[r] * rightSide[r];
+                    result[r] = _data[r] * rightSide[r];
                 }
             }
         }
@@ -696,9 +705,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 result.Clear();
 
                 // Multiply the elements in the vector with the corresponding diagonal element in this.
-                for (var r = 0; r < Data.Length; r++)
+                for (var r = 0; r < _data.Length; r++)
                 {
-                    result[r] = Data[r] * leftSide[r];
+                    result[r] = _data[r] * leftSide[r];
                 }
             }
         }
@@ -714,7 +723,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 throw new ArgumentException(Resources.ArgumentMatrixSquare);
             }
 
-            return Data.Aggregate(Complex32.One, (current, t) => current * t);
+            return _data.Aggregate(Complex32.One, (current, t) => current * t);
         }
 
         /// <summary>
@@ -727,7 +736,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         {
             // TODO: Should we return reference to array? In current implementation we return copy of array, so changes in DenseVector will
             // not influence onto diagonal elements
-            return new DenseVector((Complex32[])Data.Clone());
+            return new DenseVector((Complex32[])_data.Clone());
         }
 
         /// <summary>
@@ -795,24 +804,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public override void CopyTo(Matrix<Complex32> target)
         {
             var diagonalTarget = target as DiagonalMatrix;
-
-            if (diagonalTarget == null)
+            if (diagonalTarget != null)
             {
-                base.CopyTo(target);
+                _storage.CopyTo(diagonalTarget.Storage);
                 return;
             }
 
-            if (ReferenceEquals(this, target))
+            var denseTarget = target as DenseMatrix;
+            if (denseTarget != null)
             {
+                _storage.CopyTo(denseTarget.Storage);
                 return;
             }
 
-            if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, target, "target");
-            }
-
-            Array.Copy(Data, diagonalTarget.Data, Data.Length);
+            base.CopyTo(target);
         }
 
         /// <summary>
@@ -822,7 +827,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public override Matrix<Complex32> Transpose()
         {
             var ret = new DiagonalMatrix(ColumnCount, RowCount);
-            Array.Copy(Data, ret.Data, Data.Length);
+            Array.Copy(_data, ret._data, _data.Length);
             return ret;
         }
 
@@ -876,9 +881,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
             // Clear the result and copy the diagonal entry.
             result.Clear();
-            if (columnIndex >= rowIndex && columnIndex < rowIndex + length && columnIndex < Data.Length)
+            if (columnIndex >= rowIndex && columnIndex < rowIndex + length && columnIndex < _data.Length)
             {
-                result[columnIndex - rowIndex] = Data[columnIndex];
+                result[columnIndex - rowIndex] = _data[columnIndex];
             }
         }
 
@@ -932,9 +937,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
 
             // Clear the result and copy the diagonal entry.
             result.Clear();
-            if (rowIndex >= columnIndex && rowIndex < columnIndex + length && rowIndex < Data.Length)
+            if (rowIndex >= columnIndex && rowIndex < columnIndex + length && rowIndex < _data.Length)
             {
-                result[rowIndex - columnIndex] = Data[rowIndex];
+                result[rowIndex - columnIndex] = _data[rowIndex];
             }
         }
 
@@ -942,21 +947,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// <returns>The L1 norm of the matrix.</returns>
         public override Complex32 L1Norm()
         {
-            return Data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
+            return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
         }
 
         /// <summary>Calculates the L2 norm.</summary>
         /// <returns>The L2 norm of the matrix.</returns>   
         public override Complex32 L2Norm()
         {
-            return Data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
+            return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, t.Magnitude));
         }
 
         /// <summary>Calculates the Frobenius norm of this matrix.</summary>
         /// <returns>The Frobenius norm of this matrix.</returns>
         public override Complex32 FrobeniusNorm()
         {
-            var norm = Data.Sum(t => t.Magnitude * t.Magnitude);
+            var norm = _data.Sum(t => t.Magnitude * t.Magnitude);
             return Convert.ToSingle(Math.Sqrt(norm));
         }
 
@@ -973,7 +978,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         {
             var maxSv = float.NegativeInfinity;
             var minSv = float.PositiveInfinity;
-            foreach (var t in Data)
+            foreach (var t in _data)
             {
                 maxSv = Math.Max(maxSv, t.Magnitude);
                 minSv = Math.Min(minSv, t.Magnitude);
@@ -994,11 +999,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
 
             var inverse = (DiagonalMatrix)Clone();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                if (Data[i] != 0.0f)
+                if (_data[i] != 0.0f)
                 {
-                    inverse.Data[i] = 1.0f / Data[i];
+                    inverse._data[i] = 1.0f / _data[i];
                 }
                 else
                 {
@@ -1042,9 +1047,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
 
             result.Clear();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
         }
 
@@ -1107,9 +1112,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             }
 
             result.Clear();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
         }
 
@@ -1203,7 +1208,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 int end = Math.Min(columnCount, rowCount + columnInit);
                 for (var i = 0; columnInit + i < end; i++)
                 {
-                    result[i, columnInit + i] = Data[rowIndex + i];
+                    result[i, columnInit + i] = _data[rowIndex + i];
                 }
             }
             else if (rowIndex < columnIndex && rowIndex + rowCount > columnIndex)
@@ -1212,14 +1217,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
                 int end = Math.Min(columnCount + rowInit, rowCount);
                 for (var i = 0; rowInit + i < end; i++)
                 {
-                    result[rowInit + i, i] = Data[columnIndex + i];
+                    result[rowInit + i, i] = _data[columnIndex + i];
                 }
             }
             else
             {
                 for (var i = 0; i < Math.Min(columnCount, rowCount); i++)
                 {
-                    result[i, i] = Data[rowIndex + i];
+                    result[i, i] = _data[rowIndex + i];
                 }
             }
 
@@ -1233,9 +1238,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         public override Complex32[,] ToArray()
         {
             var result = new Complex32[RowCount, ColumnCount];
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             return result;
@@ -1384,9 +1389,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1452,9 +1457,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1516,9 +1521,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1573,7 +1578,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             var m = new DiagonalMatrix(order);
             for (var i = 0; i < order; i++)
             {
-                m.Data[i] = 1.0f;
+                m._data[i] = 1.0f;
             }
 
             return m;

src/Numerics/LinearAlgebra/Double/DenseMatrix.cs

@@ -162,6 +162,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
+        internal DenseColumnMajorMatrixStorage<double> Storage
+        {
+            get { return _storage; }
+        }
+
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>
@@ -217,7 +222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             var denseTarget = target as DenseMatrix;
             if (denseTarget != null)
             {
-                _storage.CopyTo(denseTarget._storage);
+                _storage.CopyTo(denseTarget.Storage);
                 return;
             }
 
@@ -361,7 +366,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </summary>
         public override void Clear()
         {
-            Array.Clear(_data, 0, _data.Length);
+            _storage.Clear();
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs

@@ -30,6 +30,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
     using System.Linq;
     using Generic;
     using Properties;
+    using Storage;
     using Threading;
 
     /// <summary>
@@ -44,6 +45,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
     [Serializable]
     public class DiagonalMatrix : Matrix
     {
+        readonly SparseDiagonalMatrixStorage<double> _storage;
+
+        /// <summary>
+        /// Gets the matrix's data.
+        /// </summary>
+        /// <value>The matrix's data.</value>
+        readonly double[] _data;
+
          /// <summary>
         /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class. This matrix is square with a given size.
         /// </summary>
@@ -53,7 +62,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </exception>
         public DiagonalMatrix(int order) : base(order)
         {
-            Data = new double[order];
+            _storage = new SparseDiagonalMatrixStorage<double>(order, order);
+            _data = _storage.Data;
         }
 
         /// <summary>
@@ -67,11 +77,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </param>
         public DiagonalMatrix(int rows, int columns) : base(rows, columns)
         {
-            Data = new double[Math.Min(rows, columns)];
+            _storage = new SparseDiagonalMatrixStorage<double>(rows, columns);
+            _data = _storage.Data;
         }
 
         /// <summary>
-        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class with all entries set to a particular value.
+        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class with all diagonal entries set to a particular value.
         /// </summary>
         /// <param name="rows">
         /// The number of rows.
@@ -79,13 +90,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <param name="columns">
         /// The number of columns.
         /// </param>
-        /// <param name="value">The value which we assign to each element of the matrix.</param>
+        /// <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)
         {
-            Data = new double[Math.Min(rows, columns)];
-            for (var i = 0; i < Data.Length; i++)
+            _storage = new SparseDiagonalMatrixStorage<double>(rows, columns);
+            _data = _storage.Data;
+
+            for (var i = 0; i < _data.Length; i++)
             {
-                Data[i] = value;
+                _data[i] = value;
             }
         }
 
@@ -98,7 +111,8 @@ 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)
         {
-            Data = diagonalArray;
+            _storage = new SparseDiagonalMatrixStorage<double>(rows, columns, diagonalArray);
+            _data = _storage.Data;
         }
 
         /// <summary>
@@ -110,16 +124,16 @@ 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))
         {
-            var rows = array.GetLength(0);
-            var columns = array.GetLength(1);
+            _storage = new SparseDiagonalMatrixStorage<double>(array.GetLength(0), array.GetLength(1));
+            _data = _storage.Data;
 
-            for (var i = 0; i < rows; i++)
+            for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j < columns; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     if (i == j)
                     {
-                        Data[i] = array[i, j];
+                        _data[i] = array[i, j];
                     }
                     else if (array[i, j] != 0.0 && !Double.IsNaN(array[i, j]))
                     {
@@ -129,14 +143,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
         }
 
-        /// <summary>
-        /// Gets the matrix's data.
-        /// </summary>
-        /// <value>The matrix's data.</value>
-        internal double[] Data
+        internal SparseDiagonalMatrixStorage<double> Storage
         {
-            get;
-            private set;
+            get { return _storage; }
         }
 
         /// <summary>
@@ -155,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 row == column ? _data[row] : 0.0;
         }
 
         /// <summary>
@@ -177,7 +186,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             if (row == column)
             {
-                Data[row] = value;
+                _data[row] = value;
             }
             else if (value != 0.0 && !Double.IsNaN(value))
             {
@@ -219,7 +228,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </summary>
         public override void Clear()
         {
-            Array.Clear(Data, 0, Data.Length);
+            _storage.Clear();
         }
 
         /// <summary>
@@ -246,13 +255,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return true;
             }
 
-            if (diagonalMatrix.Data.Length != Data.Length)
+            if (diagonalMatrix._data.Length != _data.Length)
             {
                 return false;
             }
 
             // If all else fails, perform element wise comparison.
-            return !Data.Where((t, i) => t != diagonalMatrix.Data[i]).Any();
+            return !_data.Where((t, i) => t != diagonalMatrix._data[i]).Any();
         }
 
         /// <summary>
@@ -263,14 +272,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// </returns>
         public override int GetHashCode()
         {
-            var hashNum = Math.Min(Data.Length, 25);
+            var hashNum = Math.Min(_data.Length, 25);
             long hash = 0;
             for (var i = 0; i < hashNum; i++)
             {
 #if PORTABLE
-                hash ^= Precision.DoubleToInt64Bits(Data[i]);
+                hash ^= Precision.DoubleToInt64Bits(_data[i]);
 #else
-                hash ^= BitConverter.DoubleToInt64Bits(Data[i]);
+                hash ^= BitConverter.DoubleToInt64Bits(_data[i]);
 #endif
             }
 
@@ -350,7 +359,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
             else
             {
-                Control.LinearAlgebraProvider.AddArrays(Data, diagOther.Data, diagResult.Data);    
+                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);    
             }
         }
 
@@ -425,7 +434,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
             else
             {
-                Control.LinearAlgebraProvider.SubtractArrays(Data, diagOther.Data, diagResult.Data);
+                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
             }
         }
 
@@ -446,12 +455,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentNullException("source");
             }
 
-            if (source.Length != Data.Length)
+            if (source.Length != _data.Length)
             {
                 throw new ArgumentException(Resources.ArgumentArraysSameLength, "source");
             }
-        
-            Buffer.BlockCopy(source, 0, Data, 0, source.Length * Constants.SizeOfDouble);
+
+            Buffer.BlockCopy(source, 0, _data, 0, source.Length * Constants.SizeOfDouble);
         }
 
         /// <summary>
@@ -473,12 +482,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 return;
             }
 
-            if (Data.Length != denseSource.Data.Length)
+            if (_data.Length != denseSource.Data.Length)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "source");
             }
 
-            Buffer.BlockCopy(denseSource.Data, 0, Data, 0, denseSource.Data.Length  * Constants.SizeOfDouble);
+            Buffer.BlockCopy(denseSource.Data, 0, _data, 0, denseSource.Data.Length * Constants.SizeOfDouble);
         }
 
         /// <summary>
@@ -514,7 +523,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                     CopyTo(diagResult);
                 }
 
-                Control.LinearAlgebraProvider.ScaleArray(scalar, Data, diagResult.Data);
+                Control.LinearAlgebraProvider.ScaleArray(scalar, _data, diagResult._data);
             }
         }
 
@@ -553,12 +562,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
             else
             {
-                var thisDataCopy = new double[r.Data.Length];
-                var otherDataCopy = new double[r.Data.Length];
-                Buffer.BlockCopy(Data, 0, thisDataCopy, 0, (r.Data.Length > Data.Length) ? Data.Length * Constants.SizeOfDouble : r.Data.Length * Constants.SizeOfDouble);
-                Buffer.BlockCopy(m.Data, 0, otherDataCopy, 0, (r.Data.Length > m.Data.Length) ? m.Data.Length * Constants.SizeOfDouble : r.Data.Length * Constants.SizeOfDouble);
+                var thisDataCopy = new double[r._data.Length];
+                var otherDataCopy = new double[r._data.Length];
+                Buffer.BlockCopy(_data, 0, thisDataCopy, 0, (r._data.Length > _data.Length) ? _data.Length * Constants.SizeOfDouble : r._data.Length * Constants.SizeOfDouble);
+                Buffer.BlockCopy(m._data, 0, otherDataCopy, 0, (r._data.Length > m._data.Length) ? m._data.Length * Constants.SizeOfDouble : r._data.Length * Constants.SizeOfDouble);
 
-                Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, r.Data);
+                Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, r._data);
             }
         }
 
@@ -635,9 +644,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 result.Clear();
 
                 // Multiply the elements in the vector with the corresponding diagonal element in this.
-                for (var r = 0; r < Data.Length; r++)
+                for (var r = 0; r < _data.Length; r++)
                 {
-                    result[r] = Data[r] * rightSide[r];
+                    result[r] = _data[r] * rightSide[r];
                 }
             }
         }
@@ -685,9 +694,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 result.Clear();
 
                 // Multiply the elements in the vector with the corresponding diagonal element in this.
-                for (var r = 0; r < Data.Length; r++)
+                for (var r = 0; r < _data.Length; r++)
                 {
-                    result[r] = Data[r] * leftSide[r];
+                    result[r] = _data[r] * leftSide[r];
                 }
             }
         }
@@ -703,7 +712,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 throw new ArgumentException(Resources.ArgumentMatrixSquare);
             }
 
-            return Data.Aggregate(1.0, (current, t) => current * t);
+            return _data.Aggregate(1.0, (current, t) => current * t);
         }
 
         /// <summary>
@@ -716,7 +725,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             // TODO: Should we return reference to array? In current implementation we return copy of array, so changes in DenseVector will
             // not influence onto diagonal elements
-            return new DenseVector((double[])Data.Clone());
+            return new DenseVector((double[])_data.Clone());
         }
 
         /// <summary>
@@ -784,24 +793,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public override void CopyTo(Matrix<double> target)
         {
             var diagonalTarget = target as DiagonalMatrix;
-
-            if (diagonalTarget == null)
+            if (diagonalTarget != null)
             {
-                base.CopyTo(target);
+                _storage.CopyTo(diagonalTarget.Storage);
                 return;
             }
 
-            if (ReferenceEquals(this, target))
+            var denseTarget = target as DenseMatrix;
+            if (denseTarget != null)
             {
+                _storage.CopyTo(denseTarget.Storage);
                 return;
             }
 
-            if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, target, "target");
-            }
-
-            Buffer.BlockCopy(Data, 0, diagonalTarget.Data, 0, Data.Length * Constants.SizeOfDouble);
+            base.CopyTo(target);
         }
 
         /// <summary>
@@ -811,7 +816,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public override Matrix<double> Transpose()
         {
             var ret = new DiagonalMatrix(ColumnCount, RowCount);
-            Buffer.BlockCopy(Data, 0, ret.Data, 0, Data.Length * Constants.SizeOfDouble);
+            Buffer.BlockCopy(_data, 0, ret._data, 0, _data.Length * Constants.SizeOfDouble);
             return ret;
         }
 
@@ -865,9 +870,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             // Clear the result and copy the diagonal entry.
             result.Clear();
-            if (columnIndex >= rowIndex && columnIndex < rowIndex + length && columnIndex < Data.Length)
+            if (columnIndex >= rowIndex && columnIndex < rowIndex + length && columnIndex < _data.Length)
             {
-                result[columnIndex - rowIndex] = Data[columnIndex];
+                result[columnIndex - rowIndex] = _data[columnIndex];
             }
         }
 
@@ -921,9 +926,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
 
             // Clear the result and copy the diagonal entry.
             result.Clear();
-            if (rowIndex >= columnIndex && rowIndex < columnIndex + length && rowIndex < Data.Length)
+            if (rowIndex >= columnIndex && rowIndex < columnIndex + length && rowIndex < _data.Length)
             {
-                result[rowIndex - columnIndex] = Data[rowIndex];
+                result[rowIndex - columnIndex] = _data[rowIndex];
             }
         }
 
@@ -931,21 +936,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// <returns>The L1 norm of the matrix.</returns>
         public override double L1Norm()
         {
-            return Data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
+            return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
         }
 
         /// <summary>Calculates the L2 norm.</summary>
         /// <returns>The L2 norm of the matrix.</returns>   
         public override double L2Norm()
         {
-            return Data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
+            return _data.Aggregate(double.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
         }
 
         /// <summary>Calculates the Frobenius norm of this matrix.</summary>
         /// <returns>The Frobenius norm of this matrix.</returns>
         public override double FrobeniusNorm()
         {
-            var norm = Data.Sum(t => t * t);
+            var norm = _data.Sum(t => t * t);
             return Math.Sqrt(norm);
         }
 
@@ -962,7 +967,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         {
             var maxSv = double.NegativeInfinity;
             var minSv = double.PositiveInfinity;
-            foreach (var t in Data)
+            foreach (var t in _data)
             {
                 maxSv = Math.Max(maxSv, Math.Abs(t));
                 minSv = Math.Min(minSv, Math.Abs(t));
@@ -983,11 +988,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             var inverse = (DiagonalMatrix)Clone();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                if (Data[i] != 0.0)
+                if (_data[i] != 0.0)
                 {
-                    inverse.Data[i] = 1.0 / Data[i];
+                    inverse._data[i] = 1.0 / _data[i];
                 }
                 else
                 {
@@ -1031,9 +1036,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             result.Clear();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
         }
 
@@ -1096,9 +1101,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             }
 
             result.Clear();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
         }
 
@@ -1192,7 +1197,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 int end = Math.Min(columnCount, rowCount + columnInit);
                 for (var i = 0; columnInit + i < end; i++)
                 {
-                    result[i, columnInit + i] = Data[rowIndex + i];
+                    result[i, columnInit + i] = _data[rowIndex + i];
                 }
             }
             else if (rowIndex < columnIndex && rowIndex + rowCount > columnIndex)
@@ -1201,14 +1206,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                 int end = Math.Min(columnCount + rowInit, rowCount);
                 for (var i = 0; rowInit + i < end; i++)
                 {
-                    result[rowInit + i, i] = Data[columnIndex + i];
+                    result[rowInit + i, i] = _data[columnIndex + i];
                 }
             }
             else
             {
                 for (var i = 0; i < Math.Min(columnCount, rowCount); i++)
                 {
-                    result[i, i] = Data[rowIndex + i];
+                    result[i, i] = _data[rowIndex + i];
                 }
             }
 
@@ -1222,9 +1227,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         public override double[,] ToArray()
         {
             var result = new double[RowCount, ColumnCount];
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             return result;
@@ -1373,9 +1378,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1441,9 +1446,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1505,9 +1510,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1567,9 +1572,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
                     CopyTo(result);
                 }
 
-                for (var index = 0; index < Data.Length; index++)
+                for (var index = 0; index < _data.Length; index++)
                 {
-                    denseResult.Data[index] %= divisor;
+                    denseResult._data[index] %= divisor;
                 }
             }
         }
@@ -1589,7 +1594,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             var m = new DiagonalMatrix(order);
             for (var i = 0; i < order; i++)
             {
-                m.Data[i] = 1.0;
+                m._data[i] = 1.0;
             }
 
             return m;

src/Numerics/LinearAlgebra/Generic/Matrix.cs

@@ -94,7 +94,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// Gets the number of columns.
         /// </summary>
         /// <value>The number of columns.</value>
-        public virtual int ColumnCount
+        public int ColumnCount
         {
             get;
             private set;
@@ -104,7 +104,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic
         /// Gets the number of rows.
         /// </summary>
         /// <value>The number of rows.</value>
-        public virtual int RowCount
+        public int RowCount
         {
             get;
             private set;

src/Numerics/LinearAlgebra/Single/DenseMatrix.cs

@@ -162,6 +162,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
+        internal DenseColumnMajorMatrixStorage<float> Storage
+        {
+            get { return _storage; }
+        }
+
         /// <summary>
         /// Gets the matrix's data.
         /// </summary>
@@ -217,7 +222,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             var denseTarget = target as DenseMatrix;
             if (denseTarget != null)
             {
-                _storage.CopyTo(denseTarget._storage);
+                _storage.CopyTo(denseTarget.Storage);
                 return;
             }
 
@@ -361,7 +366,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </summary>
         public override void Clear()
         {
-            Array.Clear(_data, 0, _data.Length);
+            _storage.Clear();
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs

@@ -30,6 +30,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
     using System.Linq;
     using Generic;
     using Properties;
+    using Storage;
     using Threading;
 
     /// <summary>
@@ -44,6 +45,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single
     [Serializable]
     public class DiagonalMatrix : Matrix
     {
+        readonly SparseDiagonalMatrixStorage<float> _storage;
+
+        /// <summary>
+        /// Gets the matrix's data.
+        /// </summary>
+        /// <value>The matrix's data.</value>
+        readonly float[] _data;
+
         /// <summary>
         /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class. This matrix is square with a given size.
         /// </summary>
@@ -54,7 +63,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public DiagonalMatrix(int order)
             : base(order)
         {
-            Data = new float[order];
+            _storage = new SparseDiagonalMatrixStorage<float>(order, order);
+            _data = _storage.Data;
         }
 
         /// <summary>
@@ -69,11 +79,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public DiagonalMatrix(int rows, int columns)
             : base(rows, columns)
         {
-            Data = new float[Math.Min(rows, columns)];
+            _storage = new SparseDiagonalMatrixStorage<float>(rows, columns);
+            _data = _storage.Data;
         }
 
         /// <summary>
-        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class with all entries set to a particular value.
+        /// Initializes a new instance of the <see cref="DiagonalMatrix"/> class with all diagonal entries set to a particular value.
         /// </summary>
         /// <param name="rows">
         /// The number of rows.
@@ -81,14 +92,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <param name="columns">
         /// The number of columns.
         /// </param>
-        /// <param name="value">The value which we assign to each element of the matrix.</param>
+        /// <param name="value">The value which we assign to each diagonal element of the matrix.</param>
         public DiagonalMatrix(int rows, int columns, float value)
             : base(rows, columns)
         {
-            Data = new float[Math.Min(rows, columns)];
-            for (var i = 0; i < Data.Length; i++)
+            _storage = new SparseDiagonalMatrixStorage<float>(rows, columns);
+            _data = _storage.Data;
+
+            for (var i = 0; i < _data.Length; i++)
             {
-                Data[i] = value;
+                _data[i] = value;
             }
         }
 
@@ -102,7 +115,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public DiagonalMatrix(int rows, int columns, float[] diagonalArray)
             : base(rows, columns)
         {
-            Data = diagonalArray;
+            _storage = new SparseDiagonalMatrixStorage<float>(rows, columns, diagonalArray);
+            _data = _storage.Data;
         }
 
         /// <summary>
@@ -115,16 +129,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public DiagonalMatrix(float[,] array)
             : this(array.GetLength(0), array.GetLength(1))
         {
-            var rows = array.GetLength(0);
-            var columns = array.GetLength(1);
+            _storage = new SparseDiagonalMatrixStorage<float>(array.GetLength(0), array.GetLength(1));
+            _data = _storage.Data;
 
-            for (var i = 0; i < rows; i++)
+            for (var i = 0; i < RowCount; i++)
             {
-                for (var j = 0; j < columns; j++)
+                for (var j = 0; j < ColumnCount; j++)
                 {
                     if (i == j)
                     {
-                        Data[i] = array[i, j];
+                        _data[i] = array[i, j];
                     }
                     else if (array[i, j] != 0.0 && !float.IsNaN(array[i, j]))
                     {
@@ -134,14 +148,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
         }
 
-        /// <summary>
-        /// Gets the matrix's data.
-        /// </summary>
-        /// <value>The matrix's data.</value>
-        internal float[] Data
+        internal SparseDiagonalMatrixStorage<float> Storage
         {
-            get;
-            private set;
+            get { return _storage; }
         }
 
         /// <summary>
@@ -160,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 row == column ? _data[row] : 0.0f;
         }
 
         /// <summary>
@@ -182,7 +191,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         {
             if (row == column)
             {
-                Data[row] = value;
+                _data[row] = value;
             }
             else if (value != 0.0 && !float.IsNaN(value))
             {
@@ -224,7 +233,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </summary>
         public override void Clear()
         {
-            Array.Clear(Data, 0, Data.Length);
+            _storage.Clear();
         }
 
         /// <summary>
@@ -251,13 +260,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 return true;
             }
 
-            if (diagonalMatrix.Data.Length != Data.Length)
+            if (diagonalMatrix._data.Length != _data.Length)
             {
                 return false;
             }
 
             // If all else fails, perform element wise comparison.
-            return !Data.Where((t, i) => t != diagonalMatrix.Data[i]).Any();
+            return !_data.Where((t, i) => t != diagonalMatrix._data[i]).Any();
         }
 
         /// <summary>
@@ -268,14 +277,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// </returns>
         public override int GetHashCode()
         {
-            var hashNum = Math.Min(Data.Length, 25);
+            var hashNum = Math.Min(_data.Length, 25);
             long hash = 0;
             for (var i = 0; i < hashNum; i++)
             {
 #if PORTABLE
-                hash ^= Precision.DoubleToInt64Bits(Data[i]);
+                hash ^= Precision.DoubleToInt64Bits(_data[i]);
 #else
-                hash ^= BitConverter.DoubleToInt64Bits(Data[i]);
+                hash ^= BitConverter.DoubleToInt64Bits(_data[i]);
 #endif
             }
 
@@ -355,7 +364,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
             else
             {
-                Control.LinearAlgebraProvider.AddArrays(Data, diagOther.Data, diagResult.Data);
+                Control.LinearAlgebraProvider.AddArrays(_data, diagOther._data, diagResult._data);
             }
         }
 
@@ -430,7 +439,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
             else
             {
-                Control.LinearAlgebraProvider.SubtractArrays(Data, diagOther.Data, diagResult.Data);
+                Control.LinearAlgebraProvider.SubtractArrays(_data, diagOther._data, diagResult._data);
             }
         }
 
@@ -451,12 +460,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 throw new ArgumentNullException("source");
             }
 
-            if (source.Length != Data.Length)
+            if (source.Length != _data.Length)
             {
                 throw new ArgumentException(Resources.ArgumentArraysSameLength, "source");
             }
 
-            Buffer.BlockCopy(source, 0, Data, 0, source.Length * Constants.SizeOfFloat);
+            Buffer.BlockCopy(source, 0, _data, 0, source.Length * Constants.SizeOfFloat);
         }
 
         /// <summary>
@@ -478,12 +487,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 return;
             }
 
-            if (Data.Length != denseSource.Data.Length)
+            if (_data.Length != denseSource.Data.Length)
             {
                 throw new ArgumentException(Resources.ArgumentVectorsSameLength, "source");
             }
 
-            Buffer.BlockCopy(denseSource.Data, 0, Data, 0, denseSource.Data.Length * Constants.SizeOfFloat);
+            Buffer.BlockCopy(denseSource.Data, 0, _data, 0, denseSource.Data.Length * Constants.SizeOfFloat);
         }
 
         /// <summary>
@@ -519,7 +528,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                     CopyTo(diagResult);
                 }
 
-                Control.LinearAlgebraProvider.ScaleArray(scalar, Data, diagResult.Data);
+                Control.LinearAlgebraProvider.ScaleArray(scalar, _data, diagResult._data);
             }
         }
 
@@ -558,12 +567,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
             else
             {
-                var thisDataCopy = new float[r.Data.Length];
-                var otherDataCopy = new float[r.Data.Length];
-                Buffer.BlockCopy(Data, 0, thisDataCopy, 0, (r.Data.Length > Data.Length) ? Data.Length * Constants.SizeOfFloat : r.Data.Length * Constants.SizeOfFloat);
-                Buffer.BlockCopy(m.Data, 0, otherDataCopy, 0, (r.Data.Length > m.Data.Length) ? m.Data.Length * Constants.SizeOfFloat : r.Data.Length * Constants.SizeOfFloat);
+                var thisDataCopy = new float[r._data.Length];
+                var otherDataCopy = new float[r._data.Length];
+                Buffer.BlockCopy(_data, 0, thisDataCopy, 0, (r._data.Length > _data.Length) ? _data.Length * Constants.SizeOfFloat : r._data.Length * Constants.SizeOfFloat);
+                Buffer.BlockCopy(m._data, 0, otherDataCopy, 0, (r._data.Length > m._data.Length) ? m._data.Length * Constants.SizeOfFloat : r._data.Length * Constants.SizeOfFloat);
 
-                Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, r.Data);
+                Control.LinearAlgebraProvider.PointWiseMultiplyArrays(thisDataCopy, otherDataCopy, r._data);
             }
         }
 
@@ -640,9 +649,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 result.Clear();
 
                 // Multiply the elements in the vector with the corresponding diagonal element in this.
-                for (var r = 0; r < Data.Length; r++)
+                for (var r = 0; r < _data.Length; r++)
                 {
-                    result[r] = Data[r] * rightSide[r];
+                    result[r] = _data[r] * rightSide[r];
                 }
             }
         }
@@ -690,9 +699,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 result.Clear();
 
                 // Multiply the elements in the vector with the corresponding diagonal element in this.
-                for (var r = 0; r < Data.Length; r++)
+                for (var r = 0; r < _data.Length; r++)
                 {
-                    result[r] = Data[r] * leftSide[r];
+                    result[r] = _data[r] * leftSide[r];
                 }
             }
         }
@@ -708,7 +717,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 throw new ArgumentException(Resources.ArgumentMatrixSquare);
             }
 
-            return Data.Aggregate(1.0f, (current, t) => current * t);
+            return _data.Aggregate(1.0f, (current, t) => current * t);
         }
 
         /// <summary>
@@ -721,7 +730,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         {
             // TODO: Should we return reference to array? In current implementation we return copy of array, so changes in DenseVector will
             // not influence onto diagonal elements
-            return new DenseVector((float[])Data.Clone());
+            return new DenseVector((float[])_data.Clone());
         }
 
         /// <summary>
@@ -789,24 +798,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public override void CopyTo(Matrix<float> target)
         {
             var diagonalTarget = target as DiagonalMatrix;
-
-            if (diagonalTarget == null)
+            if (diagonalTarget != null)
             {
-                base.CopyTo(target);
+                _storage.CopyTo(diagonalTarget.Storage);
                 return;
             }
 
-            if (ReferenceEquals(this, target))
+            var denseTarget = target as DenseMatrix;
+            if (denseTarget != null)
             {
+                _storage.CopyTo(denseTarget.Storage);
                 return;
             }
 
-            if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
-            {
-                throw DimensionsDontMatch<ArgumentException>(this, target, "target");
-            }
-
-            Buffer.BlockCopy(Data, 0, diagonalTarget.Data, 0, Data.Length * Constants.SizeOfFloat);
+            base.CopyTo(target);
         }
 
         /// <summary>
@@ -816,7 +821,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public override Matrix<float> Transpose()
         {
             var ret = new DiagonalMatrix(ColumnCount, RowCount);
-            Buffer.BlockCopy(Data, 0, ret.Data, 0, Data.Length * Constants.SizeOfFloat);
+            Buffer.BlockCopy(_data, 0, ret._data, 0, _data.Length * Constants.SizeOfFloat);
             return ret;
         }
 
@@ -870,9 +875,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
             // Clear the result and copy the diagonal entry.
             result.Clear();
-            if (columnIndex >= rowIndex && columnIndex < rowIndex + length && columnIndex < Data.Length)
+            if (columnIndex >= rowIndex && columnIndex < rowIndex + length && columnIndex < _data.Length)
             {
-                result[columnIndex - rowIndex] = Data[columnIndex];
+                result[columnIndex - rowIndex] = _data[columnIndex];
             }
         }
 
@@ -926,9 +931,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
 
             // Clear the result and copy the diagonal entry.
             result.Clear();
-            if (rowIndex >= columnIndex && rowIndex < columnIndex + length && rowIndex < Data.Length)
+            if (rowIndex >= columnIndex && rowIndex < columnIndex + length && rowIndex < _data.Length)
             {
-                result[rowIndex - columnIndex] = Data[rowIndex];
+                result[rowIndex - columnIndex] = _data[rowIndex];
             }
         }
 
@@ -936,21 +941,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// <returns>The L1 norm of the matrix.</returns>
         public override float L1Norm()
         {
-            return Data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
+            return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
         }
 
         /// <summary>Calculates the L2 norm.</summary>
         /// <returns>The L2 norm of the matrix.</returns>   
         public override float L2Norm()
         {
-            return Data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
+            return _data.Aggregate(float.NegativeInfinity, (current, t) => Math.Max(current, Math.Abs(t)));
         }
 
         /// <summary>Calculates the Frobenius norm of this matrix.</summary>
         /// <returns>The Frobenius norm of this matrix.</returns>
         public override float FrobeniusNorm()
         {
-            var norm = Data.Sum(t => t * t);
+            var norm = _data.Sum(t => t * t);
             return Convert.ToSingle(Math.Sqrt(norm));
         }
 
@@ -967,7 +972,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         {
             var maxSv = float.NegativeInfinity;
             var minSv = float.PositiveInfinity;
-            foreach (var t in Data)
+            foreach (var t in _data)
             {
                 maxSv = Math.Max(maxSv, Math.Abs(t));
                 minSv = Math.Min(minSv, Math.Abs(t));
@@ -988,11 +993,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             var inverse = (DiagonalMatrix)Clone();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                if (Data[i] != 0.0)
+                if (_data[i] != 0.0)
                 {
-                    inverse.Data[i] = 1.0f / Data[i];
+                    inverse._data[i] = 1.0f / _data[i];
                 }
                 else
                 {
@@ -1036,9 +1041,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             result.Clear();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
         }
 
@@ -1101,9 +1106,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             }
 
             result.Clear();
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
         }
 
@@ -1197,7 +1202,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 int end = Math.Min(columnCount, rowCount + columnInit);
                 for (var i = 0; columnInit + i < end; i++)
                 {
-                    result[i, columnInit + i] = Data[rowIndex + i];
+                    result[i, columnInit + i] = _data[rowIndex + i];
                 }
             }
             else if (rowIndex < columnIndex && rowIndex + rowCount > columnIndex)
@@ -1206,14 +1211,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                 int end = Math.Min(columnCount + rowInit, rowCount);
                 for (var i = 0; rowInit + i < end; i++)
                 {
-                    result[rowInit + i, i] = Data[columnIndex + i];
+                    result[rowInit + i, i] = _data[columnIndex + i];
                 }
             }
             else
             {
                 for (var i = 0; i < Math.Min(columnCount, rowCount); i++)
                 {
-                    result[i, i] = Data[rowIndex + i];
+                    result[i, i] = _data[rowIndex + i];
                 }
             }
 
@@ -1227,9 +1232,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         public override float[,] ToArray()
         {
             var result = new float[RowCount, ColumnCount];
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             return result;
@@ -1378,9 +1383,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1446,9 +1451,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1510,9 +1515,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             result.Clear();
 
             // Copy the diagonal part into the result matrix.
-            for (var i = 0; i < Data.Length; i++)
+            for (var i = 0; i < _data.Length; i++)
             {
-                result[i, i] = Data[i];
+                result[i, i] = _data[i];
             }
 
             // Copy the lower matrix into the result matrix.
@@ -1572,9 +1577,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
                     CopyTo(result);
                 }
 
-                for (var index = 0; index < Data.Length; index++)
+                for (var index = 0; index < _data.Length; index++)
                 {
-                    denseResult.Data[index] %= divisor;
+                    denseResult._data[index] %= divisor;
                 }
             }
         }
@@ -1594,7 +1599,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             var m = new DiagonalMatrix(order);
             for (var i = 0; i < order; i++)
             {
-                m.Data[i] = 1.0f;
+                m._data[i] = 1.0f;
             }
 
             return m;

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -24,6 +24,11 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             Data = data;
         }
 
+        public void Clear()
+        {
+            Array.Clear(Data, 0, Data.Length);
+        }
+
         public void CopyTo(DenseColumnMajorMatrixStorage<T> target)
         {
             if (ReferenceEquals(this, target))

src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs

@@ -0,0 +1,75 @@
+using System;
+using MathNet.Numerics.Properties;
+
+namespace MathNet.Numerics.LinearAlgebra.Storage
+{
+    internal class SparseDiagonalMatrixStorage<T>
+        where T : struct, IEquatable<T>, IFormattable
+    {
+        public int RowCount { get; private set; }
+        public int ColumnCount { get; private set; }
+        public T[] Data { get; private set; }
+
+        internal SparseDiagonalMatrixStorage(int rows, int columns)
+        {
+            RowCount = rows;
+            ColumnCount = columns;
+            Data = new T[Math.Min(rows, columns)];
+        }
+
+        internal SparseDiagonalMatrixStorage(int rows, int columns, T[] data)
+        {
+            RowCount = rows;
+            ColumnCount = columns;
+            Data = data;
+        }
+
+        public void Clear()
+        {
+            Array.Clear(Data, 0, Data.Length);
+        }
+
+        public void CopyTo(SparseDiagonalMatrixStorage<T> target)
+        {
+            if (ReferenceEquals(this, target))
+            {
+                return;
+            }
+
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
+            {
+                var message = string.Format(Resources.ArgumentMatrixDimensions2, RowCount + "x" + ColumnCount, target.RowCount + "x" + target.ColumnCount);
+                throw new ArgumentException(message, "target");
+            }
+
+            //Buffer.BlockCopy(Data, 0, target.Data, 0, Data.Length * System.Runtime.InteropServices.Marshal.SizeOf(typeof(T)));
+            Array.Copy(Data, 0, target.Data, 0, Data.Length);
+        }
+
+        public void CopyTo(DenseColumnMajorMatrixStorage<T> target)
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            if (RowCount != target.RowCount || ColumnCount != target.ColumnCount)
+            {
+                var message = string.Format(Resources.ArgumentMatrixDimensions2, RowCount + "x" + ColumnCount, target.RowCount + "x" + target.ColumnCount);
+                throw new ArgumentException(message, "target");
+            }
+
+            target.Clear();
+
+            for (int i = 0; i < Data.Length; i++)
+            {
+                target.Data[i*(target.RowCount + 1)] = Data[i];
+            }
+        }
+    }
+}

src/Numerics/Numerics.csproj

@@ -331,6 +331,7 @@
     <Compile Include="LinearAlgebra\Generic\Vector.cs" />
     <Compile Include="LinearAlgebra\Single\Vector.cs" />
     <Compile Include="LinearAlgebra\Storage\DenseColumnMajorMatrixStorage.cs" />
+    <Compile Include="LinearAlgebra\Storage\SparseDiagonalMatrixStorage.cs" />
     <Compile Include="Permutation.cs" />
     <Compile Include="Distributions\Continuous\Beta.cs" />
     <Compile Include="Distributions\Continuous\ContinuousUniform.cs" />

src/Portable/Portable.csproj

@@ -891,6 +891,9 @@
     <Compile Include="..\Numerics\LinearAlgebra\Storage\DenseColumnMajorMatrixStorage.cs">
       <Link>LinearAlgebra\Storage\DenseColumnMajorMatrixStorage.cs</Link>
     </Compile>
+    <Compile Include="..\Numerics\LinearAlgebra\Storage\SparseDiagonalMatrixStorage.cs">
+      <Link>LinearAlgebra\Storage\SparseDiagonalMatrixStorage.cs</Link>
+    </Compile>
     <Compile Include="..\Numerics\NumberTheory\IntegerTheory.cs">
       <Link>NumberTheory\IntegerTheory.cs</Link>
     </Compile>