From 5264d178d8d696379c162c27caea681972687c7a Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Tue, 28 Aug 2012 12:35:59 +0200 Subject: [PATCH] LA: sparse implementation of Kroenecker product --- .../LinearAlgebra/Complex/SparseMatrix.cs | 49 ++++++++++++++++--- .../LinearAlgebra/Complex32/SparseMatrix.cs | 49 ++++++++++++++++--- .../LinearAlgebra/Double/SparseMatrix.cs | 49 ++++++++++++++++--- .../LinearAlgebra/Single/SparseMatrix.cs | 49 ++++++++++++++++--- 4 files changed, 172 insertions(+), 24 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs index 26731193..49bf71f5 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs @@ -1057,7 +1057,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var values = _storage.Values; var valueCount = _storage.ValueCount; - for (var i = 0; i < other.RowCount; i++) + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row var startIndex = rowPointers[i]; @@ -1065,7 +1065,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex for (var j = startIndex; j < endIndex; j++) { - var resVal = values[j] * other.At(i, columnIndices[j]); + var resVal = values[j]*other.At(i, columnIndices[j]); if (!resVal.IsZero()) { result.At(i, columnIndices[j], resVal); @@ -1088,7 +1088,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex var values = _storage.Values; var valueCount = _storage.ValueCount; - for (var i = 0; i < other.RowCount; i++) + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row var startIndex = rowPointers[i]; @@ -1096,10 +1096,47 @@ namespace MathNet.Numerics.LinearAlgebra.Complex for (var j = startIndex; j < endIndex; j++) { - var resVal = values[j] / other.At(i, columnIndices[j]); - if (!resVal.IsZero()) + if (!values[j].IsZero()) { - result.At(i, columnIndices[j], resVal); + result.At(i, columnIndices[j], values[j]/other.At(i, columnIndices[j])); + } + } + } + } + + public override void KroneckerProduct(Matrix other, Matrix result) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != (RowCount*other.RowCount) || result.ColumnCount != (ColumnCount*other.ColumnCount)) + { + throw DimensionsDontMatch(this, other, result); + } + + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + for (var i = 0; i < RowCount; i++) + { + // Get the begin / end index for the current row + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; + + for (var j = startIndex; j < endIndex; j++) + { + if (!values[j].IsZero()) + { + result.SetSubMatrix(i*other.RowCount, other.RowCount, columnIndices[j]*other.ColumnCount, other.ColumnCount, values[j]*other); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs index ed3b25b6..1c46f3ac 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs @@ -1056,7 +1056,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var values = _storage.Values; var valueCount = _storage.ValueCount; - for (var i = 0; i < other.RowCount; i++) + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row var startIndex = rowPointers[i]; @@ -1064,7 +1064,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 for (var j = startIndex; j < endIndex; j++) { - var resVal = values[j] * other.At(i, columnIndices[j]); + var resVal = values[j]*other.At(i, columnIndices[j]); if (!resVal.IsZero()) { result.At(i, columnIndices[j], resVal); @@ -1087,7 +1087,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 var values = _storage.Values; var valueCount = _storage.ValueCount; - for (var i = 0; i < other.RowCount; i++) + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row var startIndex = rowPointers[i]; @@ -1095,10 +1095,47 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 for (var j = startIndex; j < endIndex; j++) { - var resVal = values[j] / other.At(i, columnIndices[j]); - if (!resVal.IsZero()) + if (!values[j].IsZero()) { - result.At(i, columnIndices[j], resVal); + result.At(i, columnIndices[j], values[j]/other.At(i, columnIndices[j])); + } + } + } + } + + public override void KroneckerProduct(Matrix other, Matrix result) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != (RowCount*other.RowCount) || result.ColumnCount != (ColumnCount*other.ColumnCount)) + { + throw DimensionsDontMatch(this, other, result); + } + + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + for (var i = 0; i < RowCount; i++) + { + // Get the begin / end index for the current row + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; + + for (var j = startIndex; j < endIndex; j++) + { + if (!values[j].IsZero()) + { + result.SetSubMatrix(i*other.RowCount, other.RowCount, columnIndices[j]*other.ColumnCount, other.ColumnCount, values[j]*other); } } } diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs index 47d6cd96..24b83055 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs @@ -1055,7 +1055,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double var values = _storage.Values; var valueCount = _storage.ValueCount; - for (var i = 0; i < other.RowCount; i++) + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row var startIndex = rowPointers[i]; @@ -1063,7 +1063,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var j = startIndex; j < endIndex; j++) { - var resVal = values[j] * other.At(i, columnIndices[j]); + var resVal = values[j]*other.At(i, columnIndices[j]); if (resVal != 0d) { result.At(i, columnIndices[j], resVal); @@ -1086,7 +1086,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double var values = _storage.Values; var valueCount = _storage.ValueCount; - for (var i = 0; i < other.RowCount; i++) + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row var startIndex = rowPointers[i]; @@ -1094,10 +1094,47 @@ namespace MathNet.Numerics.LinearAlgebra.Double for (var j = startIndex; j < endIndex; j++) { - var resVal = values[j] / other.At(i, columnIndices[j]); - if (resVal != 0d) + if (values[j] != 0d) { - result.At(i, columnIndices[j], resVal); + result.At(i, columnIndices[j], values[j]/other.At(i, columnIndices[j])); + } + } + } + } + + public override void KroneckerProduct(Matrix other, Matrix result) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != (RowCount*other.RowCount) || result.ColumnCount != (ColumnCount*other.ColumnCount)) + { + throw DimensionsDontMatch(this, other, result); + } + + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + for (var i = 0; i < RowCount; i++) + { + // Get the begin / end index for the current row + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; + + for (var j = startIndex; j < endIndex; j++) + { + if (values[j] != 0d) + { + result.SetSubMatrix(i*other.RowCount, other.RowCount, columnIndices[j]*other.ColumnCount, other.ColumnCount, values[j]*other); } } } diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs index 514a796b..271ec7a6 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs @@ -1054,7 +1054,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single var values = _storage.Values; var valueCount = _storage.ValueCount; - for (var i = 0; i < other.RowCount; i++) + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row var startIndex = rowPointers[i]; @@ -1062,7 +1062,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single for (var j = startIndex; j < endIndex; j++) { - var resVal = values[j] * other.At(i, columnIndices[j]); + var resVal = values[j]*other.At(i, columnIndices[j]); if (resVal != 0f) { result.At(i, columnIndices[j], resVal); @@ -1085,7 +1085,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single var values = _storage.Values; var valueCount = _storage.ValueCount; - for (var i = 0; i < other.RowCount; i++) + for (var i = 0; i < RowCount; i++) { // Get the begin / end index for the current row var startIndex = rowPointers[i]; @@ -1093,10 +1093,47 @@ namespace MathNet.Numerics.LinearAlgebra.Single for (var j = startIndex; j < endIndex; j++) { - var resVal = values[j] / other.At(i, columnIndices[j]); - if (resVal != 0f) + if (values[j] != 0f) { - result.At(i, columnIndices[j], resVal); + result.At(i, columnIndices[j], values[j]/other.At(i, columnIndices[j])); + } + } + } + } + + public override void KroneckerProduct(Matrix other, Matrix result) + { + if (other == null) + { + throw new ArgumentNullException("other"); + } + + if (result == null) + { + throw new ArgumentNullException("result"); + } + + if (result.RowCount != (RowCount*other.RowCount) || result.ColumnCount != (ColumnCount*other.ColumnCount)) + { + throw DimensionsDontMatch(this, other, result); + } + + var rowPointers = _storage.RowPointers; + var columnIndices = _storage.ColumnIndices; + var values = _storage.Values; + var valueCount = _storage.ValueCount; + + for (var i = 0; i < RowCount; i++) + { + // Get the begin / end index for the current row + var startIndex = rowPointers[i]; + var endIndex = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount; + + for (var j = startIndex; j < endIndex; j++) + { + if (values[j] != 0f) + { + result.SetSubMatrix(i*other.RowCount, other.RowCount, columnIndices[j]*other.ColumnCount, other.ColumnCount, values[j]*other); } } }