LA: transpose at storage level, more efficient sparse implementation (via wo80)

src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs

@@ -428,42 +428,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return Control.LinearAlgebraProvider.MatrixNorm(Norm.FrobeniusNorm, _rowCount, _columnCount, _values);
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<Complex> Transpose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i];
-                }
-            }
-            return ret;
-        }
-
-        /// <summary>
-        /// Returns the conjugate transpose of this matrix.
-        /// </summary>
-        /// <returns>The conjugate transpose of this matrix.</returns>
-        public override Matrix<Complex> ConjugateTranspose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i].Conjugate();
-                }
-            }
-            return ret;
-        }
-
         /// <summary>
         /// Negate each element of this matrix and place the results into the result matrix.
         /// </summary>

src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs

@@ -728,17 +728,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return new DenseVector(_data).Clone();
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<Complex> Transpose()
-        {
-            var ret = new DiagonalMatrix(ColumnCount, RowCount);
-            Array.Copy(_data, ret._data, _data.Length);
-            return ret;
-        }
-
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
         public override double L1Norm()

src/Numerics/LinearAlgebra/Complex/Matrix.cs

@@ -110,16 +110,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// Returns the conjugate transpose of this matrix.
         /// </summary>
         /// <returns>The conjugate transpose of this matrix.</returns>
-        public override Matrix<Complex> ConjugateTranspose()
+        public override sealed Matrix<Complex> ConjugateTranspose()
         {
-            var ret = Build.SameAs(this, ColumnCount, RowCount);
+            var ret = Transpose();
-            for (var j = 0; j < ColumnCount; j++)
+            ret.MapInplace(c => c.Conjugate(), forceMapZeros: false);
-            {
-                for (var i = 0; i < RowCount; i++)
-                {
-                    ret.At(j, i, At(i, j).Conjugate());
-                }
-            }
             return ret;
         }
 

src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@@ -648,42 +648,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             DoMultiply(-1, result);
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<Complex> Transpose()
-        {
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-
-            var ret = new SparseCompressedRowMatrixStorage<Complex>(ColumnCount, RowCount)
-                {
-                    ColumnIndices = new int[_storage.ValueCount],
-                    Values = new Complex[_storage.ValueCount]
-                };
-
-            // Do an 'inverse' CopyTo iterate over the rows
-            for (var i = 0; i < RowCount; i++)
-            {
-                var startIndex = rowPointers[i];
-                var endIndex = rowPointers[i + 1];
-
-                if (startIndex == endIndex)
-                {
-                    continue;
-                }
-
-                for (var j = startIndex; j < endIndex; j++)
-                {
-                    ret.At(columnIndices[j], i, values[j]);
-                }
-            }
-
-            return new SparseMatrix(ret);
-        }
-
         /// <summary>Calculates the induced infinity norm of this matrix.</summary>
         /// <returns>The maximum absolute row sum of the matrix.</returns>
         public override double InfinityNorm()

src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs

@@ -455,42 +455,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             base.DoConjugate(result);
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<Complex32> Transpose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i];
-                }
-            }
-            return ret;
-        }
-
-        /// <summary>
-        /// Returns the conjugate transpose of this matrix.
-        /// </summary>
-        /// <returns>The conjugate transpose of this matrix.</returns>
-        public override Matrix<Complex32> ConjugateTranspose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i].Conjugate();
-                }
-            }
-            return ret;
-        }
-
         /// <summary>
         /// Add a scalar to each element of the matrix and stores the result in the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs

@@ -722,17 +722,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return new DenseVector(_data).Clone();
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<Complex32> Transpose()
-        {
-            var ret = new DiagonalMatrix(ColumnCount, RowCount);
-            Array.Copy(_data, ret._data, _data.Length);
-            return ret;
-        }
-
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
         public override double L1Norm()

src/Numerics/LinearAlgebra/Complex32/Matrix.cs

@@ -104,16 +104,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// Returns the conjugate transpose of this matrix.
         /// </summary>
         /// <returns>The conjugate transpose of this matrix.</returns>
-        public override Matrix<Complex32> ConjugateTranspose()
+        public override sealed Matrix<Complex32> ConjugateTranspose()
         {
-            var ret = Build.SameAs(this, ColumnCount, RowCount);
+            var ret = Transpose();
-            for (var j = 0; j < ColumnCount; j++)
+            ret.MapInplace(c => c.Conjugate(), forceMapZeros: false);
-            {
-                for (var i = 0; i < RowCount; i++)
-                {
-                    ret.At(j, i, At(i, j).Conjugate());
-                }
-            }
             return ret;
         }
 

src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@@ -643,42 +643,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             DoMultiply(-1, result);
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<Complex32> Transpose()
-        {
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-
-            var ret = new SparseCompressedRowMatrixStorage<Complex32>(ColumnCount, RowCount)
-                {
-                    ColumnIndices = new int[_storage.ValueCount],
-                    Values = new Complex32[_storage.ValueCount]
-                };
-
-            // Do an 'inverse' CopyTo iterate over the rows
-            for (var i = 0; i < RowCount; i++)
-            {
-                var startIndex = rowPointers[i];
-                var endIndex = rowPointers[i + 1];
-
-                if (startIndex == endIndex)
-                {
-                    continue;
-                }
-
-                for (var j = startIndex; j < endIndex; j++)
-                {
-                    ret.At(columnIndices[j], i, values[j]);
-                }
-            }
-
-            return new SparseMatrix(ret);
-        }
-
         /// <summary>Calculates the induced infinity norm of this matrix.</summary>
         /// <returns>The maximum absolute row sum of the matrix.</returns>
         public override double InfinityNorm()

src/Numerics/LinearAlgebra/Double/DenseMatrix.cs

@@ -436,25 +436,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             base.DoNegate(result);
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<double> Transpose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i];
-                }
-            }
-
-            return ret;
-        }
-
         /// <summary>
         /// Add a scalar to each element of the matrix and stores the result in the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs

@@ -572,17 +572,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return new DenseVector(_data).Clone();
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<double> Transpose()
-        {
-            var ret = new DiagonalMatrix(ColumnCount, RowCount);
-            Buffer.BlockCopy(_data, 0, ret._data, 0, _data.Length * Constants.SizeOfDouble);
-            return ret;
-        }
-
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
         public override double L1Norm()

src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@@ -641,43 +641,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             DoMultiply(-1, result);
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<double> Transpose()
-        {
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-
-            var ret = new SparseCompressedRowMatrixStorage<double>(ColumnCount, RowCount)
-                {
-                    ColumnIndices = new int[_storage.ValueCount],
-                    Values = new double[_storage.ValueCount]
-                };
-
-            // Do an 'inverse' CopyTo iterate over the rows
-            for (var i = 0; i < RowCount; i++)
-            {
-                var startIndex = rowPointers[i];
-                var endIndex = rowPointers[i + 1];
-
-                if (startIndex == endIndex)
-                {
-                    // Begin and end are equal. There are no values in the row, Move to the next row
-                    continue;
-                }
-
-                for (var j = startIndex; j < endIndex; j++)
-                {
-                    ret.At(columnIndices[j], i, values[j]);
-                }
-            }
-
-            return new SparseMatrix(ret);
-        }
-
         /// <summary>Calculates the induced infinity norm of this matrix.</summary>
         /// <returns>The maximum absolute row sum of the matrix.</returns>
         public override double InfinityNorm()

src/Numerics/LinearAlgebra/Matrix.cs

@@ -980,16 +980,10 @@ namespace MathNet.Numerics.LinearAlgebra
         /// Returns the transpose of this matrix.
         /// </summary>
         /// <returns>The transpose of this matrix.</returns>
-        public virtual Matrix<T> Transpose()
+        public Matrix<T> Transpose()
         {
             var result = Build.SameAs(this, ColumnCount, RowCount);
-            for (var j = 0; j < ColumnCount; j++)
+            Storage.TransposeToUnchecked(result.Storage, skipClearing:true);
-            {
-                for (var i = 0; i < RowCount; i++)
-                {
-                    result.At(j, i, At(i, j));
-                }
-            }
             return result;
         }
 

src/Numerics/LinearAlgebra/Single/DenseMatrix.cs

@@ -436,25 +436,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             base.DoNegate(result);
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<float> Transpose()
-        {
-            var ret = new DenseMatrix(_columnCount, _rowCount);
-            for (var j = 0; j < _columnCount; j++)
-            {
-                var index = j * _rowCount;
-                for (var i = 0; i < _rowCount; i++)
-                {
-                    ret._values[(i * _columnCount) + j] = _values[index + i];
-                }
-            }
-
-            return ret;
-        }
-
         /// <summary>
         /// Add a scalar to each element of the matrix and stores the result in the result vector.
         /// </summary>

src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs

@@ -572,17 +572,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return new DenseVector(_data).Clone();
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<float> Transpose()
-        {
-            var ret = new DiagonalMatrix(ColumnCount, RowCount);
-            Buffer.BlockCopy(_data, 0, ret._data, 0, _data.Length * Constants.SizeOfFloat);
-            return ret;
-        }
-
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
         public override double L1Norm()

src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@@ -641,45 +641,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             DoMultiply(-1, result);
         }
 
-        /// <summary>
-        /// Returns the transpose of this matrix.
-        /// </summary>
-        /// <returns>The transpose of this matrix.</returns>
-        public override Matrix<float> Transpose()
-        {
-            var rowPointers = _storage.RowPointers;
-            var columnIndices = _storage.ColumnIndices;
-            var values = _storage.Values;
-
-            var ret = new SparseCompressedRowMatrixStorage<float>(ColumnCount, RowCount)
-                {
-                    ColumnIndices = new int[_storage.ValueCount],
-                    Values = new float[_storage.ValueCount]
-                };
-
-            // Do an 'inverse' CopyTo iterate over the rows
-            for (var i = 0; i < RowCount; i++)
-            {
-                // Get the begin / end index for the current row
-                var startIndex = rowPointers[i];
-                var endIndex = rowPointers[i + 1];
-
-                // Get the values for the current row
-                if (startIndex == endIndex)
-                {
-                    // Begin and end are equal. There are no values in the row, Move to the next row
-                    continue;
-                }
-
-                for (var j = startIndex; j < endIndex; j++)
-                {
-                    ret.At(columnIndices[j], i, values[j]);
-                }
-            }
-
-            return new SparseMatrix(ret);
-        }
-
         /// <summary>Calculates the induced infinity norm of this matrix.</summary>
         /// <returns>The maximum absolute row sum of the matrix.</returns>
         public override double InfinityNorm()

src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@@ -432,6 +432,72 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
+        // TRANSPOSE
+
+        internal override void TransposeToUnchecked(MatrixStorage<T> target, bool skipClearing = false)
+        {
+            var denseTarget = target as DenseColumnMajorMatrixStorage<T>;
+            if (denseTarget != null)
+            {
+                TransposeToUnchecked(denseTarget);
+                return;
+            }
+
+            var sparseTarget = target as SparseCompressedRowMatrixStorage<T>;
+            if (sparseTarget != null)
+            {
+                TransposeToUnchecked(sparseTarget);
+                return;
+            }
+
+            // FALL BACK
+
+            for (int j = 0, offset = 0; j < ColumnCount; j++, offset += RowCount)
+            {
+                for (int i = 0; i < RowCount; i++)
+                {
+                    target.At(j, i, Data[i + offset]);
+                }
+            }
+        }
+
+        void TransposeToUnchecked(DenseColumnMajorMatrixStorage<T> target)
+        {
+            for (var j = 0; j < ColumnCount; j++)
+            {
+                var index = j * RowCount;
+                for (var i = 0; i < RowCount; i++)
+                {
+                    target.Data[(i * ColumnCount) + j] = Data[index + i];
+                }
+            }
+        }
+
+        void TransposeToUnchecked(SparseCompressedRowMatrixStorage<T> target)
+        {
+            var rowPointers = target.RowPointers;
+            var columnIndices = new List<int>();
+            var values = new List<T>();
+
+            for (int j = 0; j < ColumnCount; j++)
+            {
+                rowPointers[j] = values.Count;
+                var index = j * RowCount;
+                for (int i = 0; i < RowCount; i++)
+                {
+                    if (!Zero.Equals(Data[index + i]))
+                    {
+                        values.Add(Data[index + i]);
+                        columnIndices.Add(i);
+                    }
+                }
+            }
+
+            rowPointers[ColumnCount] = values.Count;
+            target.ColumnIndices = columnIndices.ToArray();
+            target.Values = values.ToArray();
+        }
+
         // EXTRACT
 
         public override T[] ToRowMajorArray()

src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs

@@ -494,6 +494,13 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
+        // TRANSPOSE
+
+        internal override void TransposeToUnchecked(MatrixStorage<T> target, bool skipClearing = false)
+        {
+            CopyToUnchecked(target, skipClearing);
+        }
+
         // EXTRACT
 
         public override T[] ToRowMajorArray()

src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs

@@ -384,6 +384,40 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
+        // TRANSPOSE
+
+        public void TransposeTo(MatrixStorage<T> target, bool skipClearing = false)
+        {
+            if (target == null)
+            {
+                throw new ArgumentNullException("target");
+            }
+
+            if (ReferenceEquals(this, target))
+            {
+                throw new NotSupportedException("In-place transpose is not supported.");
+            }
+
+            if (RowCount != target.ColumnCount || ColumnCount != target.RowCount)
+            {
+                var message = string.Format(Resources.ArgumentMatrixDimensions2, RowCount + "x" + ColumnCount, target.RowCount + "x" + target.ColumnCount);
+                throw new ArgumentException(message, "target");
+            }
+
+            TransposeToUnchecked(target, skipClearing);
+        }
+
+        internal virtual void TransposeToUnchecked(MatrixStorage<T> target, bool skipClearing = false)
+        {
+            for (int j = 0; j < ColumnCount; j++)
+            {
+                for (int i = 0; i < RowCount; i++)
+                {
+                    target.At(j, i, At(i, j));
+                }
+            }
+        }
+
         // EXTRACT
 
         public virtual T[] ToRowMajorArray()

src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@@ -851,6 +851,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
                 target.Clear();
             }
 
+            // TODO: proper implementation
+
             if (ValueCount != 0)
             {
                 for (int row = 0; row < RowCount; row++)
@@ -1008,6 +1010,105 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
             }
         }
 
+        // TRANSPOSE
+
+        internal override void TransposeToUnchecked(MatrixStorage<T> target, bool skipClearing = false)
+        {
+            var sparseTarget = target as SparseCompressedRowMatrixStorage<T>;
+            if (sparseTarget != null)
+            {
+                TransposeToUnchecked(sparseTarget);
+                return;
+            }
+
+            var denseTarget = target as DenseColumnMajorMatrixStorage<T>;
+            if (denseTarget != null)
+            {
+                TransposeToUnchecked(denseTarget, skipClearing);
+                return;
+            }
+
+            // FALL BACK
+
+            if (!skipClearing)
+            {
+                target.Clear();
+            }
+
+            if (ValueCount != 0)
+            {
+                for (int row = 0; row < RowCount; row++)
+                {
+                    var startIndex = RowPointers[row];
+                    var endIndex = RowPointers[row + 1];
+                    for (var j = startIndex; j < endIndex; j++)
+                    {
+                        target.At(ColumnIndices[j], row, Values[j]);
+                    }
+                }
+            }
+        }
+
+        void TransposeToUnchecked(SparseCompressedRowMatrixStorage<T> target)
+        {
+            target.Values = new T[ValueCount];
+            target.ColumnIndices = new int[ValueCount];
+            var cx = target.Values;
+            var cp = target.RowPointers;
+            var ci = target.ColumnIndices;
+
+            // Column counts
+            int[] w = new int[ColumnCount];
+            for (int p = 0; p < RowPointers[RowCount]; p++)
+            {
+                w[ColumnIndices[p]]++;
+            }
+
+            // Column pointers
+            int nz = 0;
+            for (int i = 0; i < ColumnCount; i++)
+            {
+                cp[i] = nz;
+                nz += w[i];
+                w[i] = cp[i];
+            }
+            cp[ColumnCount] = nz;
+
+            for (int i = 0; i < RowCount; i++)
+            {
+                for (int p = RowPointers[i]; p < RowPointers[i + 1]; p++)
+                {
+                    int j = w[ColumnIndices[p]]++;
+
+                    // Place A(i,j) as entry C(j,i)
+                    ci[j] = i;
+                    cx[j] = Values[p];
+                }
+            }
+        }
+
+        void TransposeToUnchecked(DenseColumnMajorMatrixStorage<T> target, bool skipClearing)
+        {
+            if (!skipClearing)
+            {
+                target.Clear();
+            }
+
+            if (ValueCount != 0)
+            {
+                for (int row = 0; row < RowCount; row++)
+                {
+                    var targetIndex = row * ColumnCount;
+                    var startIndex = RowPointers[row];
+                    var endIndex = RowPointers[row + 1];
+                    for (var j = startIndex; j < endIndex; j++)
+                    {
+                        target.Data[targetIndex + ColumnIndices[j]] = Values[j];
+                    }
+                }
+            }
+        }
+
         // EXTRACT
 
         public override T[] ToRowMajorArray()