LA Storage: diagonal sub matrix

src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs

@@ -1174,68 +1174,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// is not positive.</exception>
         public override Matrix<Complex> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
         {
-            if (rowIndex >= RowCount || rowIndex < 0)
+            // TODO: if rowIndex == columnIndex, use a diagonal matrix instead of a sparse one
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "rowCount");
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "columnCount");
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
 
-            if (colMax > ColumnCount)
+            var storage = new SparseCompressedRowMatrixStorage<Complex>(rowCount, columnCount, Complex.Zero);
-            {
+            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, true);
-                throw new ArgumentOutOfRangeException("columnCount");
+            return new SparseMatrix(storage);
-            }
-
-            var result = new SparseMatrix(rowCount, columnCount);
-
-            if (rowIndex > columnIndex && columnIndex + columnCount > rowIndex)
-            {
-                int columnInit = rowIndex - columnIndex;
-                int end = Math.Min(columnCount, rowCount + columnInit);
-                for (var i = 0; columnInit + i < end; i++)
-                {
-                    result[i, columnInit + i] = _data[rowIndex + i];
-                }
-            }
-            else if (rowIndex < columnIndex && rowIndex + rowCount > columnIndex)
-            {
-                int rowInit = columnIndex - rowIndex;
-                int end = Math.Min(columnCount + rowInit, rowCount);
-                for (var i = 0; rowInit + i < end; 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];
-                }
-            }
-
-            return result;
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@@ -56,6 +56,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             get { return _storage.ValueCount; }
         }
 
+        internal SparseMatrix(SparseCompressedRowMatrixStorage<Complex> storage)
+            : base(storage.RowCount, storage.ColumnCount)
+        {
+            _storage = storage;
+        }
+
         /// <summary>
         /// Initializes a new instance of the <see cref="SparseMatrix"/> class.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs

@@ -1179,68 +1179,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// is not positive.</exception>
         public override Matrix<Complex32> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
         {
-            if (rowIndex >= RowCount || rowIndex < 0)
+            // TODO: if rowIndex == columnIndex, use a diagonal matrix instead of a sparse one
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "rowCount");
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "columnCount");
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
 
-            if (colMax > ColumnCount)
+            var storage = new SparseCompressedRowMatrixStorage<Complex32>(rowCount, columnCount, Complex32.Zero);
-            {
+            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, true);
-                throw new ArgumentOutOfRangeException("columnCount");
+            return new SparseMatrix(storage);
-            }
-
-            var result = new SparseMatrix(rowCount, columnCount);
-
-            if (rowIndex > columnIndex && columnIndex + columnCount > rowIndex)
-            {
-                int columnInit = rowIndex - columnIndex;
-                int end = Math.Min(columnCount, rowCount + columnInit);
-                for (var i = 0; columnInit + i < end; i++)
-                {
-                    result[i, columnInit + i] = _data[rowIndex + i];
-                }
-            }
-            else if (rowIndex < columnIndex && rowIndex + rowCount > columnIndex)
-            {
-                int rowInit = columnIndex - rowIndex;
-                int end = Math.Min(columnCount + rowInit, rowCount);
-                for (var i = 0; rowInit + i < end; 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];
-                }
-            }
-
-            return result;
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@@ -56,6 +56,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             get { return _storage.ValueCount; }
         }
 
+        internal SparseMatrix(SparseCompressedRowMatrixStorage<Complex32> storage)
+            : base(storage.RowCount, storage.ColumnCount)
+        {
+            _storage = storage;
+        }
+
         /// <summary>
         /// Initializes a new instance of the <see cref="SparseMatrix"/> class.
         /// </summary>

src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs

@@ -1168,68 +1168,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
         /// is not positive.</exception>
         public override Matrix<double> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
         {
-            if (rowIndex >= RowCount || rowIndex < 0)
+            // TODO: if rowIndex == columnIndex, use a diagonal matrix instead of a sparse one
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "rowCount");
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "columnCount");
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
 
-            if (colMax > ColumnCount)
+            var storage = new SparseCompressedRowMatrixStorage<double>(rowCount, columnCount, 0d);
-            {
+            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, true);
-                throw new ArgumentOutOfRangeException("columnCount");
+            return new SparseMatrix(storage);
-            }
-
-            var result = new SparseMatrix(rowCount, columnCount);
-
-            if (rowIndex > columnIndex && columnIndex + columnCount > rowIndex)
-            {
-                int columnInit = rowIndex - columnIndex;
-                int end = Math.Min(columnCount, rowCount + columnInit);
-                for (var i = 0; columnInit + i < end; i++)
-                {
-                    result[i, columnInit + i] = _data[rowIndex + i];
-                }
-            }
-            else if (rowIndex < columnIndex && rowIndex + rowCount > columnIndex)
-            {
-                int rowInit = columnIndex - rowIndex;
-                int end = Math.Min(columnCount + rowInit, rowCount);
-                for (var i = 0; rowInit + i < end; 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];
-                }
-            }
-
-            return result;
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -55,6 +55,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             get { return _storage.ValueCount; }
         }
 
+        internal SparseMatrix(SparseCompressedRowMatrixStorage<double> storage)
+            : base(storage.RowCount, storage.ColumnCount)
+        {
+            _storage = storage;
+        }
+
         /// <summary>
         /// Initializes a new instance of the <see cref="SparseMatrix"/> class.
         /// </summary>

src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs

@@ -1173,68 +1173,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
         /// is not positive.</exception>
         public override Matrix<float> SubMatrix(int rowIndex, int rowCount, int columnIndex, int columnCount)
         {
-            if (rowIndex >= RowCount || rowIndex < 0)
+            // TODO: if rowIndex == columnIndex, use a diagonal matrix instead of a sparse one
-            {
-                throw new ArgumentOutOfRangeException("rowIndex");
-            }
-
-            if (columnIndex >= ColumnCount || columnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("columnIndex");
-            }
-
-            if (rowCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "rowCount");
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "columnCount");
-            }
-
-            var colMax = columnIndex + columnCount;
-            var rowMax = rowIndex + rowCount;
-
-            if (rowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
 
-            if (colMax > ColumnCount)
+            var storage = new SparseCompressedRowMatrixStorage<float>(rowCount, columnCount, 0f);
-            {
+            _storage.CopySubMatrixTo(storage, rowIndex, 0, rowCount, columnIndex, 0, columnCount, true);
-                throw new ArgumentOutOfRangeException("columnCount");
+            return new SparseMatrix(storage);
-            }
-
-            var result = new SparseMatrix(rowCount, columnCount);
-
-            if (rowIndex > columnIndex && columnIndex + columnCount > rowIndex)
-            {
-                int columnInit = rowIndex - columnIndex;
-                int end = Math.Min(columnCount, rowCount + columnInit);
-                for (var i = 0; columnInit + i < end; i++)
-                {
-                    result[i, columnInit + i] = _data[rowIndex + i];
-                }
-            }
-            else if (rowIndex < columnIndex && rowIndex + rowCount > columnIndex)
-            {
-                int rowInit = columnIndex - rowIndex;
-                int end = Math.Min(columnCount + rowInit, rowCount);
-                for (var i = 0; rowInit + i < end; 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];
-                }
-            }
-
-            return result;
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@@ -55,6 +55,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             get { return _storage.ValueCount; }
         }
 
+        internal SparseMatrix(SparseCompressedRowMatrixStorage<float> storage)
+            : base(storage.RowCount, storage.ColumnCount)
+        {
+            _storage = storage;
+        }
+
         /// <summary>
         /// Initializes a new instance of the <see cref="SparseMatrix"/> class.
         /// </summary>

src/Numerics/LinearAlgebra/Storage/ArgumentValidation.cs

@@ -0,0 +1,77 @@
+using System;
+using MathNet.Numerics.Properties;
+
+namespace MathNet.Numerics.LinearAlgebra.Storage
+{
+    // ReSharper disable UnusedParameter.Global
+    internal static class ArgumentValidation
+    {
+        public static void CopySubMatrixTo(
+            int sourceRowCount, int sourceColumnCount,
+            int targetRowCount, int targetColumnCount,
+            int sourceRowIndex, int targetRowIndex, int rowCount,
+            int sourceColumnIndex, int targetColumnIndex, int columnCount)
+        {
+            if (rowCount < 1)
+            {
+                throw new ArgumentException(Resources.ArgumentMustBePositive, "rowCount");
+            }
+
+            if (columnCount < 1)
+            {
+                throw new ArgumentException(Resources.ArgumentMustBePositive, "columnCount");
+            }
+
+            // Verify Source
+
+            if (sourceRowIndex >= sourceRowCount || sourceRowIndex < 0)
+            {
+                throw new ArgumentOutOfRangeException("sourceRowIndex");
+            }
+
+            if (sourceColumnIndex >= sourceColumnCount || sourceColumnIndex < 0)
+            {
+                throw new ArgumentOutOfRangeException("sourceColumnIndex");
+            }
+
+            var sourceRowMax = sourceRowIndex + rowCount;
+            var sourceColumnMax = sourceColumnIndex + columnCount;
+
+            if (sourceRowMax > sourceRowCount)
+            {
+                throw new ArgumentOutOfRangeException("rowCount");
+            }
+
+            if (sourceColumnMax > sourceColumnCount)
+            {
+                throw new ArgumentOutOfRangeException("columnCount");
+            }
+
+            // Verify Target
+
+            if (targetRowIndex >= targetRowCount || targetRowIndex < 0)
+            {
+                throw new ArgumentOutOfRangeException("targetRowIndex");
+            }
+
+            if (targetColumnIndex >= targetColumnCount || targetColumnIndex < 0)
+            {
+                throw new ArgumentOutOfRangeException("targetColumnIndex");
+            }
+
+            var targetRowMax = targetRowIndex + rowCount;
+            var targetColumnMax = targetColumnIndex + columnCount;
+
+            if (targetRowMax > targetRowCount)
+            {
+                throw new ArgumentOutOfRangeException("rowCount");
+            }
+
+            if (targetColumnMax > targetColumnCount)
+            {
+                throw new ArgumentOutOfRangeException("columnCount");
+            }
+        }
+    }
+    // ReSharper restore UnusedParameter.Global
+}

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -121,43 +121,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             int sourceRowIndex, int targetRowIndex, int rowCount,
             int sourceColumnIndex, int targetColumnIndex, int columnCount)
         {
-            if (rowCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "rowCount");
-            }
-
-            if (columnCount < 1)
-            {
-                throw new ArgumentException(Resources.ArgumentMustBePositive, "columnCount");
-            }
-
-            // Verify Source
-
-            if (sourceRowIndex >= RowCount || sourceRowIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("sourceRowIndex");
-            }
-
-            if (sourceColumnIndex >= ColumnCount || sourceColumnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("sourceColumnIndex");
-            }
-
-            var sourceRowMax = sourceRowIndex + rowCount;
-            var sourceColumnMax = sourceColumnIndex + columnCount;
-
-            if (sourceRowMax > RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
-
-            if (sourceColumnMax > ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("columnCount");
-            }
-
-            // Verify Target
-
             if (target == null)
             {
                 throw new ArgumentNullException("target");
@@ -168,32 +131,12 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 throw new NotSupportedException();
             }
 
-            if (targetRowIndex >= target.RowCount || targetRowIndex < 0)
+            ArgumentValidation.CopySubMatrixTo(RowCount, ColumnCount,
-            {
+                target.RowCount, target.ColumnCount,
-                throw new ArgumentOutOfRangeException("targetRowIndex");
+                sourceRowIndex, targetRowIndex, rowCount,
-            }
+                sourceColumnIndex, targetColumnIndex, columnCount);
-
-            if (targetColumnIndex >= target.ColumnCount || targetColumnIndex < 0)
-            {
-                throw new ArgumentOutOfRangeException("targetColumnIndex");
-            }
-
-            var targetRowMax = targetRowIndex + rowCount;
-            var targetColumnMax = targetColumnIndex + columnCount;
-
-            if (targetRowMax > target.RowCount)
-            {
-                throw new ArgumentOutOfRangeException("rowCount");
-            }
-
-            if (targetColumnMax > target.ColumnCount)
-            {
-                throw new ArgumentOutOfRangeException("columnCount");
-            }
-
-            // Copy
 
-            for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnMax; j++, jj++)
+            for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnIndex + columnCount; j++, jj++)
             {
                 //Buffer.BlockCopy(Data, j*RowCount + sourceRowIndex, target.Data, jj*target.RowCount + targetRowIndex, rowCount * System.Runtime.InteropServices.Marshal.SizeOf(typeof(T)));
                 Array.Copy(Data, j*RowCount + sourceRowIndex, target.Data, jj*target.RowCount + targetRowIndex, rowCount);

src/Numerics/LinearAlgebra/Storage/SparseDiagonalMatrixStorage.cs

@@ -174,5 +174,113 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 target.Data[i*(target.RowCount + 1)] = Data[i];
             }
         }
+
+        public void CopySubMatrixTo(DenseColumnMajorMatrixStorage<T> target,
+            int sourceRowIndex, int targetRowIndex, int rowCount,
+            int sourceColumnIndex, int targetColumnIndex, int columnCount,
+            bool skipClearing = false)
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            ArgumentValidation.CopySubMatrixTo(RowCount, ColumnCount,
+                target.RowCount, target.ColumnCount,
+                sourceRowIndex, targetRowIndex, rowCount,
+                sourceColumnIndex, targetColumnIndex, columnCount);
+
+            if (!skipClearing)
+            {
+                target.Clear();
+            }
+
+            if (sourceRowIndex > sourceColumnIndex && sourceColumnIndex + columnCount > sourceRowIndex)
+            {
+                // column by column, but skip resulting zero columns at the beginning
+
+                int columnInit = sourceRowIndex - sourceColumnIndex;
+                int offset = (columnInit + targetColumnIndex) * target.RowCount + targetRowIndex;
+                int step = target.RowCount + 1;
+                int end = Math.Min(columnCount - columnInit, rowCount) + sourceRowIndex;
+
+                for (int i = sourceRowIndex, j = offset; i < end; i++, j += step)
+                {
+                    target.Data[j] = Data[i];
+                }
+            }
+            else if (sourceRowIndex < sourceColumnIndex && sourceRowIndex + rowCount > sourceColumnIndex)
+            {
+                // row by row, but skip resulting zero rows at the beginning
+
+                int rowInit = sourceColumnIndex - sourceRowIndex;
+                int offset = targetColumnIndex*target.RowCount + rowInit + targetRowIndex;
+                int step = target.RowCount + 1;
+                int end = Math.Min(columnCount, rowCount - rowInit) + sourceColumnIndex;
+
+                for (int i = sourceColumnIndex, j = offset; i < end; i++, j += step)
+                {
+                    target.Data[j] = Data[i];
+                }
+            }
+            else
+            {
+                int offset = targetColumnIndex*target.RowCount + targetRowIndex;
+                int step = target.RowCount + 1;
+                var end = Math.Min(columnCount, rowCount) + sourceRowIndex;
+
+                for (int i = sourceRowIndex, j = offset; i < end; i++, j += step)
+                {
+                    target.Data[j] = Data[i];
+                }
+            }
+        }
+
+        public void CopySubMatrixTo(SparseCompressedRowMatrixStorage<T> target,
+            int sourceRowIndex, int targetRowIndex, int rowCount,
+            int sourceColumnIndex, int targetColumnIndex, int columnCount,
+            bool skipClearing = false)
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            ArgumentValidation.CopySubMatrixTo(RowCount, ColumnCount,
+                target.RowCount, target.ColumnCount,
+                sourceRowIndex, targetRowIndex, rowCount,
+                sourceColumnIndex, targetColumnIndex, columnCount);
+
+            if (!skipClearing)
+            {
+                target.Clear();
+            }
+
+            if (sourceRowIndex > sourceColumnIndex && sourceColumnIndex + columnCount > sourceRowIndex)
+            {
+                // column by column, but skip resulting zero columns at the beginning
+                int columnInit = sourceRowIndex - sourceColumnIndex;
+                for (var i = 0; i < Math.Min(columnCount - columnInit, rowCount); i++)
+                {
+                    target.At(i + targetRowIndex, columnInit + i + targetColumnIndex, Data[sourceRowIndex + i]);
+                }
+            }
+            else if (sourceRowIndex < sourceColumnIndex && sourceRowIndex + rowCount > sourceColumnIndex)
+            {
+                // row by row, but skip resulting zero rows at the beginning
+                int rowInit = sourceColumnIndex - sourceRowIndex;
+                for (var i = 0; i < Math.Min(columnCount, rowCount - rowInit); i++)
+                {
+                    target.At(rowInit + i + targetRowIndex, i + targetColumnIndex, Data[sourceColumnIndex + i]);
+                }
+            }
+            else
+            {
+                for (var i = 0; i < Math.Min(columnCount, rowCount); i++)
+                {
+                    target.At(i + targetRowIndex, i + targetColumnIndex, Data[sourceRowIndex + i]);
+                }
+            }
+        }
     }
 }

src/Numerics/Numerics.csproj

@@ -122,6 +122,7 @@
     <Compile Include="Distributions\Multivariate\InverseWishart.cs" />
     <Compile Include="Distributions\Multivariate\MatrixNormal.cs" />
     <Compile Include="Distributions\Multivariate\Wishart.cs" />
+    <Compile Include="LinearAlgebra\Storage\ArgumentValidation.cs" />
     <Compile Include="LinearAlgebra\Complex32\ExtensionMethods.cs" />
     <Compile Include="LinearAlgebra\Complex32\Factorization\Cholesky.cs" />
     <Compile Include="LinearAlgebra\Complex32\DenseMatrix.cs" />

src/Portable/Portable.csproj

@@ -888,6 +888,9 @@
     <Compile Include="..\Numerics\LinearAlgebra\Single\Vector.cs">
       <Link>LinearAlgebra\Single\Vector.cs</Link>
     </Compile>
+    <Compile Include="..\Numerics\LinearAlgebra\Storage\ArgumentValidation.cs">
+      <Link>LinearAlgebra\Storage\ArgumentValidation.cs</Link>
+    </Compile>
     <Compile Include="..\Numerics\LinearAlgebra\Storage\DenseColumnMajorMatrixStorage.cs">
       <Link>LinearAlgebra\Storage\DenseColumnMajorMatrixStorage.cs</Link>
     </Compile>