Browse Source

LA: Matrix MapSubMatrixIndexedTo

provider
Christoph Ruegg 13 years ago
parent
commit
1671a9ec32
  1. 23
      src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
  2. 17
      src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
  3. 23
      src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
  4. 17
      src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
  5. 23
      src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
  6. 17
      src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
  7. 4
      src/Numerics/LinearAlgebra/Matrix.cs
  8. 23
      src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
  9. 17
      src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
  10. 86
      src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs
  11. 295
      src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs
  12. 19
      src/Numerics/LinearAlgebra/Storage/MatrixStorage.Validation.cs
  13. 87
      src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs
  14. 392
      src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs
  15. 260
      src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Map.cs
  16. 4
      src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Reform.cs
  17. 1
      src/UnitTests/UnitTests.csproj

23
src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs

@ -435,14 +435,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return;
}
if (RowCount == ColumnCount)
if (ColumnCount == RowCount)
{
// TODO: Map is rather generic, we can do better
other.MapIndexed((i, j, x) => _data[i]*x, result, false);
return;
other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false);
}
else
{
result.Clear();
other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true);
}
base.DoMultiply(other, result);
}
/// <summary>
@ -577,7 +578,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return;
}
base.DoTransposeThisAndMultiply(other, result);
if (ColumnCount == RowCount)
{
other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false);
}
else
{
result.Clear();
other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true);
}
}
/// <summary>

17
src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs

@ -931,12 +931,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return;
}
var diagonalOther = other as DiagonalMatrix;
if (diagonalOther != null && sparseResult != null && other.RowCount == other.ColumnCount)
var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex>;
if (diagonalOther != null && sparseResult != null)
{
var diagonal = ((DiagonalMatrixStorage<Complex>)other.Storage).Data;
// TODO: Map is rather generic, we can do better
MapIndexed((i, j, x) => x*diagonal[j], result, false);
var diagonal = diagonalOther.Data;
if (other.ColumnCount == other.RowCount)
{
Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], false, false);
}
else
{
result.Storage.Clear();
Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, ColumnCount, false, true);
}
return;
}

23
src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs

@ -429,14 +429,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return;
}
if (RowCount == ColumnCount)
if (ColumnCount == RowCount)
{
// TODO: Map is rather generic, we can do better
other.MapIndexed((i, j, x) => _data[i]*x, result, false);
return;
other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false);
}
else
{
result.Clear();
other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true);
}
base.DoMultiply(other, result);
}
/// <summary>
@ -571,7 +572,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return;
}
base.DoTransposeThisAndMultiply(other, result);
if (ColumnCount == RowCount)
{
other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false);
}
else
{
result.Clear();
other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true);
}
}
/// <summary>

17
src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs

@ -925,12 +925,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return;
}
var diagonalOther = other as DiagonalMatrix;
if (diagonalOther != null && sparseResult != null && other.RowCount == other.ColumnCount)
var diagonalOther = other.Storage as DiagonalMatrixStorage<Complex32>;
if (diagonalOther != null && sparseResult != null)
{
var diagonal = ((DiagonalMatrixStorage<Complex32>)other.Storage).Data;
// TODO: Map is rather generic, we can do better
MapIndexed((i, j, x) => x*diagonal[j], result, false);
var diagonal = diagonalOther.Data;
if (other.ColumnCount == other.RowCount)
{
Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], false, false);
}
else
{
result.Storage.Clear();
Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, ColumnCount, false, true);
}
return;
}

23
src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs

@ -409,14 +409,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return;
}
if (RowCount == ColumnCount)
if (ColumnCount == RowCount)
{
// TODO: Map is rather generic, we can do better
other.MapIndexed((i, j, x) => _data[i]*x, result, false);
return;
other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false);
}
else
{
result.Clear();
other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true);
}
base.DoMultiply(other, result);
}
/// <summary>
@ -505,7 +506,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return;
}
base.DoTransposeThisAndMultiply(other, result);
if (ColumnCount == RowCount)
{
other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false);
}
else
{
result.Clear();
other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true);
}
}
/// <summary>

17
src/Numerics/LinearAlgebra/Double/SparseMatrix.cs

@ -927,12 +927,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return;
}
var diagonalOther = other as DiagonalMatrix;
if (diagonalOther != null && sparseResult != null && other.RowCount == other.ColumnCount)
var diagonalOther = other.Storage as DiagonalMatrixStorage<double>;
if (diagonalOther != null && sparseResult != null)
{
var diagonal = ((DiagonalMatrixStorage<double>)other.Storage).Data;
// TODO: Map is rather generic, we can do better
MapIndexed((i, j, x) => x*diagonal[j], result, false);
var diagonal = diagonalOther.Data;
if (other.ColumnCount == other.RowCount)
{
Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], false, false);
}
else
{
result.Storage.Clear();
Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, ColumnCount, false, true);
}
return;
}

4
src/Numerics/LinearAlgebra/Matrix.cs

@ -1518,7 +1518,7 @@ namespace MathNet.Numerics.LinearAlgebra
public Matrix<TU> Map<TU>(Func<T, TU> f, bool forceMapZeros = false)
where TU : struct, IEquatable<TU>, IFormattable
{
var result = Matrix<TU>.Build.SameAs(this);
var result = Matrix<TU>.Build.SameAs(this, RowCount, ColumnCount, fullyMutable: forceMapZeros);
Storage.MapToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true);
return result;
}
@ -1532,7 +1532,7 @@ namespace MathNet.Numerics.LinearAlgebra
public Matrix<TU> MapIndexed<TU>(Func<int, int, T, TU> f, bool forceMapZeros = false)
where TU : struct, IEquatable<TU>, IFormattable
{
var result = Matrix<TU>.Build.SameAs(this);
var result = Matrix<TU>.Build.SameAs(this, RowCount, ColumnCount, fullyMutable: forceMapZeros);
Storage.MapIndexedToUnchecked(result.Storage, f, forceMapZeros, skipClearing: true);
return result;
}

23
src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs

@ -409,14 +409,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return;
}
if (RowCount == ColumnCount)
if (ColumnCount == RowCount)
{
// TODO: Map is rather generic, we can do better
other.MapIndexed((i, j, x) => _data[i]*x, result, false);
return;
other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false);
}
else
{
result.Clear();
other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true);
}
base.DoMultiply(other, result);
}
/// <summary>
@ -505,7 +506,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return;
}
base.DoTransposeThisAndMultiply(other, result);
if (ColumnCount == RowCount)
{
other.Storage.MapIndexedTo(result.Storage, (i, j, x) => x*_data[i], false, false);
}
else
{
result.Clear();
other.Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*_data[i], 0, 0, other.RowCount, 0, 0, other.ColumnCount, false, true);
}
}
/// <summary>

17
src/Numerics/LinearAlgebra/Single/SparseMatrix.cs

@ -932,12 +932,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return;
}
var diagonalOther = other as DiagonalMatrix;
if (diagonalOther != null && sparseResult != null && other.RowCount == other.ColumnCount)
var diagonalOther = other.Storage as DiagonalMatrixStorage<float>;
if (diagonalOther != null && sparseResult != null)
{
var diagonal = ((DiagonalMatrixStorage<float>)other.Storage).Data;
// TODO: Map is rather generic, we can do better
MapIndexed((i, j, x) => x*diagonal[j], result, false);
var diagonal = diagonalOther.Data;
if (other.ColumnCount == other.RowCount)
{
Storage.MapIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], false, false);
}
else
{
result.Storage.Clear();
Storage.MapSubMatrixIndexedTo(result.Storage, (i, j, x) => x*diagonal[j], 0, 0, RowCount, 0, 0, ColumnCount, false, true);
}
return;
}

86
src/Numerics/LinearAlgebra/Storage/DenseColumnMajorMatrixStorage.cs

@ -365,9 +365,18 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
return;
}
// TODO: Proper Sparse Implementation
// FALL BACK
base.CopySubMatrixToUnchecked(target, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing);
for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnIndex + columnCount; j++, jj++)
{
int index = sourceRowIndex + j*RowCount;
for (int ii = targetRowIndex; ii < targetRowIndex + rowCount; ii++)
{
target.At(ii, jj, Data[index++]);
}
}
}
void CopySubMatrixToUnchecked(DenseColumnMajorMatrixStorage<T> target,
@ -503,6 +512,33 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// FUNCTIONAL COMBINATORS
public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
{
CommonParallel.For(0, Data.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
Data[i] = f(Data[i]);
}
});
}
public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
{
CommonParallel.For(0, ColumnCount, Math.Max(4096/RowCount, 32), (a, b) =>
{
int index = a*RowCount;
for (int j = a; j < b; j++)
{
for (int i = 0; i < RowCount; i++)
{
Data[index] = f(i, j, Data[index]);
index++;
}
}
});
}
internal override void MapToUnchecked<TU>(MatrixStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
{
var denseTarget = target as DenseColumnMajorMatrixStorage<TU>;
@ -530,17 +566,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
{
CommonParallel.For(0, Data.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
Data[i] = f(Data[i]);
}
});
}
internal override void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
{
var denseTarget = target as DenseColumnMajorMatrixStorage<TU>;
@ -573,20 +598,41 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
internal override void MapSubMatrixIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount,
bool forceMapZeros = false, bool skipClearing = false)
{
CommonParallel.For(0, ColumnCount, Math.Max(4096/RowCount, 32), (a, b) =>
var denseTarget = target as DenseColumnMajorMatrixStorage<TU>;
if (denseTarget != null)
{
int index = a*RowCount;
for (int j = a; j < b; j++)
CommonParallel.For(0, columnCount, Math.Max(4096/rowCount, 32), (a, b) =>
{
for (int i = 0; i < RowCount; i++)
for (int j = a; j < b; j++)
{
Data[index] = f(i, j, Data[index]);
index++;
int sourceIndex = sourceRowIndex + (j + sourceColumnIndex)*RowCount;
int targetIndex = targetRowIndex + (j + targetColumnIndex)*target.RowCount;
for (int i = 0; i < rowCount; i++)
{
denseTarget.Data[targetIndex++] = f(targetRowIndex + i, targetColumnIndex + j, Data[sourceIndex++]);
}
}
});
return;
}
// TODO: Proper Sparse Implementation
// FALL BACK
for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnIndex + columnCount; j++, jj++)
{
int index = sourceRowIndex + j*RowCount;
for (int ii = targetRowIndex; ii < targetRowIndex + rowCount; ii++)
{
target.At(ii, jj, f(ii, jj, Data[index++]));
}
});
}
}
}
}

295
src/Numerics/LinearAlgebra/Storage/DiagonalMatrixStorage.cs

@ -349,16 +349,40 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
return;
}
var sparseTarget = target as SparseCompressedRowMatrixStorage<T>;
if (sparseTarget != null)
{
CopySubMatrixToUnchecked(sparseTarget, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing);
return;
}
// TODO: Proper Sparse Implementation
// FALL BACK
base.CopySubMatrixToUnchecked(target, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, skipClearing);
if (!skipClearing)
{
target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
}
if (sourceRowIndex == sourceColumnIndex)
{
for (var i = 0; i < Math.Min(columnCount, rowCount); i++)
{
target.At(targetRowIndex + i, targetColumnIndex + i, Data[sourceRowIndex + i]);
}
}
else 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(targetRowIndex + i, columnInit + targetColumnIndex + i, 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 + targetRowIndex + i, targetColumnIndex + i, Data[sourceColumnIndex + i]);
}
}
}
void CopySubMatrixToUnchecked(DiagonalMatrixStorage<T> target,
@ -436,45 +460,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
void CopySubMatrixToUnchecked(SparseCompressedRowMatrixStorage<T> target,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount,
bool skipClearing)
{
if (!skipClearing)
{
target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
}
if (sourceRowIndex == sourceColumnIndex)
{
for (var i = 0; i < Math.Min(columnCount, rowCount); i++)
{
target.At(i + targetRowIndex, i + targetColumnIndex, Data[sourceRowIndex + i]);
}
}
else 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: all zero, nop
}
// ROW COPY
internal override void CopySubRowToUnchecked(VectorStorage<T> target, int rowIndex,
@ -589,6 +574,38 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// FUNCTIONAL COMBINATORS
public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
{
if (forceMapZeros)
{
throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix");
}
CommonParallel.For(0, Data.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
Data[i] = f(Data[i]);
}
});
}
public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
{
if (forceMapZeros)
{
throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix");
}
CommonParallel.For(0, Data.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
Data[i] = f(i, i, Data[i]);
}
});
}
internal override void MapToUnchecked<TU>(MatrixStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
{
var processZeros = forceMapZeros || !Zero.Equals(f(Zero));
@ -637,25 +654,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
{
if (forceMapZeros)
{
throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix");
}
CommonParallel.For(0, Data.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
Data[i] = f(Data[i]);
}
});
}
internal override void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
{
var processZeros = forceMapZeros || !Zero.Equals(f(0, 0, Zero));
var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero));
var diagonalTarget = target as DiagonalMatrixStorage<TU>;
if (diagonalTarget != null)
@ -701,20 +702,176 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
internal override void MapSubMatrixIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount,
bool forceMapZeros = false, bool skipClearing = false)
{
if (forceMapZeros)
var diagonalTarget = target as DiagonalMatrixStorage<TU>;
if (diagonalTarget != null)
{
MapSubMatrixIndexedToUnchecked(diagonalTarget, f, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, forceMapZeros);
return;
}
var denseTarget = target as DenseColumnMajorMatrixStorage<TU>;
if (denseTarget != null)
{
MapSubMatrixIndexedToUnchecked(denseTarget, f, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, forceMapZeros, skipClearing);
return;
}
// TODO: Proper Sparse Implementation
// FALL BACK
if (!skipClearing)
{
target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
}
if (sourceRowIndex == sourceColumnIndex)
{
int targetRow = targetRowIndex;
int targetColumn = targetColumnIndex;
for (var i = 0; i < Math.Min(columnCount, rowCount); i++)
{
target.At(targetRow, targetColumn, f(targetRow, targetColumn, Data[sourceRowIndex + i]));
targetRow++;
targetColumn++;
}
}
else if (sourceRowIndex > sourceColumnIndex && sourceColumnIndex + columnCount > sourceRowIndex)
{
// column by column, but skip resulting zero columns at the beginning
int columnInit = sourceRowIndex - sourceColumnIndex;
int targetRow = targetRowIndex;
int targetColumn = targetColumnIndex + columnInit;
for (var i = 0; i < Math.Min(columnCount - columnInit, rowCount); i++)
{
target.At(targetRow, targetColumn, f(targetRow, targetColumn, Data[sourceRowIndex + i]));
targetRow++;
targetColumn++;
}
}
else if (sourceRowIndex < sourceColumnIndex && sourceRowIndex + rowCount > sourceColumnIndex)
{
// row by row, but skip resulting zero rows at the beginning
int rowInit = sourceColumnIndex - sourceRowIndex;
int targetRow = targetRowIndex + rowInit;
int targetColumn = targetColumnIndex;
for (var i = 0; i < Math.Min(columnCount, rowCount - rowInit); i++)
{
target.At(targetRow, targetColumn, f(targetRow, targetColumn, Data[sourceColumnIndex + i]));
targetRow++;
targetColumn++;
}
}
}
void MapSubMatrixIndexedToUnchecked<TU>(DiagonalMatrixStorage<TU> target, Func<int, int, T, TU> f,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount,
bool forceMapZeros = false)
where TU : struct, IEquatable<TU>, IFormattable
{
var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero));
if (processZeros || sourceRowIndex - sourceColumnIndex != targetRowIndex - targetColumnIndex)
{
throw new NotSupportedException("Cannot map non-zero off-diagonal values into a diagonal matrix");
}
CommonParallel.For(0, Data.Length, 4096, (a, b) =>
var beginInclusive = Math.Max(sourceRowIndex, sourceColumnIndex);
var count = Math.Min(sourceRowIndex + rowCount, sourceColumnIndex + columnCount) - beginInclusive;
if (count > 0)
{
for (int i = a; i < b; i++)
var beginTarget = Math.Max(targetRowIndex, targetColumnIndex);
CommonParallel.For(0, count, 4096, (a, b) =>
{
Data[i] = f(i, i, Data[i]);
int targetIndex = beginTarget + a;
for (int i = a; i < b; i++)
{
target.Data[targetIndex] = f(targetIndex, targetIndex, Data[beginInclusive + i]);
targetIndex++;
}
});
}
}
void MapSubMatrixIndexedToUnchecked<TU>(DenseColumnMajorMatrixStorage<TU> target, Func<int, int, T, TU> f,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount,
bool forceMapZeros = false, bool skipClearing = false)
where TU : struct, IEquatable<TU>, IFormattable
{
var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero));
if (!skipClearing && !processZeros)
{
target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
}
if (processZeros)
{
CommonParallel.For(0, columnCount, Math.Max(4096/rowCount, 32), (a, b) =>
{
int sourceColumn = sourceColumnIndex + a;
int targetColumn = targetColumnIndex + a;
for (int j = a; j < b; j++)
{
int targetIndex = targetRowIndex + (j + targetColumnIndex)*target.RowCount;
int sourceRow = sourceRowIndex;
int targetRow = targetRowIndex;
for (int i = 0; i < rowCount; i++)
{
target.Data[targetIndex++] = f(targetRow++, targetColumn, sourceRow++ == sourceColumn ? Data[sourceColumn] : Zero);
}
sourceColumn++;
targetColumn++;
}
});
}
else
{
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 count = Math.Min(columnCount - columnInit, rowCount);
for (int k = 0, j = offset; k < count; j += step, k++)
{
target.Data[j] = f(targetRowIndex + k, targetColumnIndex + columnInit + k, Data[sourceRowIndex + k]);
}
}
});
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 count = Math.Min(columnCount, rowCount - rowInit);
for (int k = 0, j = offset; k < count; j += step, k++)
{
target.Data[j] = f(targetRowIndex + rowInit + k, targetColumnIndex + k, Data[sourceColumnIndex + k]);
}
}
else
{
int offset = targetColumnIndex*target.RowCount + targetRowIndex;
int step = target.RowCount + 1;
var count = Math.Min(columnCount, rowCount);
for (int k = 0, j = offset; k < count; j += step, k++)
{
target.Data[j] = f(targetRowIndex + k, targetColumnIndex + k, Data[sourceRowIndex + k]);
}
}
}
}
}
}

19
src/Numerics/LinearAlgebra/Storage/MatrixStorage.Validation.cs

@ -4,7 +4,7 @@
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2013 Math.NET
// Copyright (c) 2009-2014 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
@ -34,6 +34,7 @@ using MathNet.Numerics.Properties;
namespace MathNet.Numerics.LinearAlgebra.Storage
{
// ReSharper disable UnusedParameter.Local
public partial class MatrixStorage<T>
{
void ValidateRange(int row, int column)
@ -49,9 +50,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
void ValidateSubMatrixRange(MatrixStorage<T> target,
void ValidateSubMatrixRange<TU>(MatrixStorage<TU> target,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount)
where TU : struct, IEquatable<TU>, IFormattable
{
if (rowCount < 1)
{
@ -114,7 +116,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
void ValidateRowRange(VectorStorage<T> target, int rowIndex)
void ValidateRowRange<TU>(VectorStorage<TU> target, int rowIndex)
where TU : struct, IEquatable<TU>, IFormattable
{
if (rowIndex >= RowCount || rowIndex < 0)
{
@ -127,7 +130,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
void ValidateColumnRange(VectorStorage<T> target, int columnIndex)
void ValidateColumnRange<TU>(VectorStorage<TU> target, int columnIndex)
where TU : struct, IEquatable<TU>, IFormattable
{
if (columnIndex >= ColumnCount || columnIndex < 0)
{
@ -140,8 +144,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
void ValidateSubRowRange(VectorStorage<T> target, int rowIndex,
void ValidateSubRowRange<TU>(VectorStorage<TU> target, int rowIndex,
int sourceColumnIndex, int targetColumnIndex, int columnCount)
where TU : struct, IEquatable<TU>, IFormattable
{
if (columnCount < 1)
{
@ -178,8 +183,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
void ValidateSubColumnRange(VectorStorage<T> target, int columnIndex,
void ValidateSubColumnRange<TU>(VectorStorage<TU> target, int columnIndex,
int sourceRowIndex, int targetRowIndex, int rowCount)
where TU : struct, IEquatable<TU>, IFormattable
{
if (rowCount < 1)
{
@ -216,5 +222,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
}
// ReSharper restore UnusedParameter.Local
}

87
src/Numerics/LinearAlgebra/Storage/MatrixStorage.cs

@ -388,10 +388,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
public virtual T[] ToRowMajorArray()
{
var ret = new T[RowCount * ColumnCount];
var ret = new T[RowCount*ColumnCount];
for (int i = 0; i < RowCount; i++)
{
var offset = i * ColumnCount;
var offset = i*ColumnCount;
for (int j = 0; j < ColumnCount; j++)
{
ret[offset + j] = At(i, j);
@ -402,10 +402,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
public virtual T[] ToColumnMajorArray()
{
var ret = new T[RowCount * ColumnCount];
var ret = new T[RowCount*ColumnCount];
for (int j = 0; j < ColumnCount; j++)
{
var offset = j * RowCount;
var offset = j*RowCount;
for (int i = 0; i < RowCount; i++)
{
ret[offset + i] = At(i, j);
@ -416,7 +416,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
public virtual T[,] ToArray()
{
var ret = new T[RowCount,ColumnCount];
var ret = new T[RowCount, ColumnCount];
for (int i = 0; i < RowCount; i++)
{
for (int j = 0; j < ColumnCount; j++)
@ -483,6 +483,28 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// FUNCTIONAL COMBINATORS
public virtual void MapInplace(Func<T, T> f, bool forceMapZeros = false)
{
for (int i = 0; i < RowCount; i++)
{
for (int j = 0; j < ColumnCount; j++)
{
At(i, j, f(At(i, j)));
}
}
}
public virtual void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
{
for (int i = 0; i < RowCount; i++)
{
for (int j = 0; j < ColumnCount; j++)
{
At(i, j, f(i, j, At(i, j)));
}
}
}
public void MapTo<TU>(MatrixStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
where TU : struct, IEquatable<TU>, IFormattable
{
@ -512,17 +534,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
public virtual void MapInplace(Func<T, T> f, bool forceMapZeros = false)
{
for (int i = 0; i < RowCount; i++)
{
for (int j = 0; j < ColumnCount; j++)
{
At(i, j, f(At(i, j)));
}
}
}
public void MapIndexedTo<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
where TU : struct, IEquatable<TU>, IFormattable
{
@ -543,22 +554,54 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
internal virtual void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
where TU : struct, IEquatable<TU>, IFormattable
{
for (int i = 0; i < RowCount; i++)
for (int j = 0; j < ColumnCount; j++)
{
for (int j = 0; j < ColumnCount; j++)
for (int i = 0; i < RowCount; i++)
{
target.At(i, j, f(i, j, At(i, j)));
}
}
}
public virtual void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
public void MapSubMatrixIndexedTo<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount,
bool forceMapZeros = false, bool skipClearing = false)
where TU : struct, IEquatable<TU>, IFormattable
{
for (int i = 0; i < RowCount; i++)
if (target == null)
{
for (int j = 0; j < ColumnCount; j++)
throw new ArgumentNullException("target");
}
if (rowCount == 0 || columnCount == 0)
{
return;
}
if (ReferenceEquals(this, target))
{
throw new NotSupportedException();
}
ValidateSubMatrixRange(target,
sourceRowIndex, targetRowIndex, rowCount,
sourceColumnIndex, targetColumnIndex, columnCount);
MapSubMatrixIndexedToUnchecked(target, f, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, forceMapZeros, skipClearing);
}
internal virtual void MapSubMatrixIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount,
bool forceMapZeros = false, bool skipClearing = false)
where TU : struct, IEquatable<TU>, IFormattable
{
for (int j = sourceColumnIndex, jj = targetColumnIndex; j < sourceColumnIndex + columnCount; j++, jj++)
{
for (int i = sourceRowIndex, ii = targetRowIndex; i < sourceRowIndex + rowCount; i++, ii++)
{
At(i, j, f(i, j, At(i, j)));
target.At(ii, jj, f(ii, jj, At(i, j)));
}
}
}

392
src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

@ -926,7 +926,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
var columnIndices = new List<int>(ValueCount);
var rowPointers = target.RowPointers;
for (int i = sourceRowIndex, row = 0; i < sourceRowIndex + rowCount; i++, row++)
for (int i = sourceRowIndex; i < sourceRowIndex + rowCount; i++)
{
rowPointers[i + rowOffset] = values.Count;
@ -934,13 +934,13 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
var endIndex = RowPointers[i + 1];
// note: we might be able to replace this loop with Array.Copy (perf)
for (int j = startIndex; j < endIndex; j++)
for (int k = startIndex; k < endIndex; k++)
{
// check if the column index is in the range
if ((ColumnIndices[j] >= sourceColumnIndex) && (ColumnIndices[j] < sourceColumnIndex + columnCount))
if ((ColumnIndices[k] >= sourceColumnIndex) && (ColumnIndices[k] < sourceColumnIndex + columnCount))
{
values.Add(Values[j]);
columnIndices.Add(ColumnIndices[j] + columnOffset);
values.Add(Values[k]);
columnIndices.Add(ColumnIndices[k] + columnOffset);
}
}
}
@ -963,7 +963,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
// NOTE: potential for more efficient implementation
for (int i = sourceRowIndex, row = 0; i < sourceRowIndex + rowCount; i++, row++)
for (int i = sourceRowIndex, row = 0; row < rowCount; i++, row++)
{
var startIndex = RowPointers[i];
var endIndex = RowPointers[i + 1];
@ -1118,6 +1118,112 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// FUNCTIONAL COMBINATORS
public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
{
if (forceMapZeros || !Zero.Equals(f(Zero)))
{
var newRowPointers = RowPointers;
var newColumnIndices = new List<int>(ColumnIndices.Length);
var newValues = new List<T>(Values.Length);
int k = 0;
for (int row = 0; row < RowCount; row++)
{
newRowPointers[row] = newValues.Count;
for (int col = 0; col < ColumnCount; col++)
{
var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(Values[k++]) : f(Zero);
if (!Zero.Equals(item))
{
newValues.Add(item);
newColumnIndices.Add(col);
}
}
}
ColumnIndices = newColumnIndices.ToArray();
Values = newValues.ToArray();
newRowPointers[RowCount] = newValues.Count;
}
else
{
// we can safely do this in-place:
int nonZero = 0;
for (int row = 0; row < RowCount; row++)
{
var startIndex = RowPointers[row];
var endIndex = RowPointers[row + 1];
RowPointers[row] = nonZero;
for (var j = startIndex; j < endIndex; j++)
{
var item = f(Values[j]);
if (!Zero.Equals(item))
{
Values[nonZero] = item;
ColumnIndices[nonZero] = ColumnIndices[j];
nonZero++;
}
}
}
Array.Resize(ref ColumnIndices, nonZero);
Array.Resize(ref Values, nonZero);
RowPointers[RowCount] = nonZero;
}
}
public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
{
if (forceMapZeros || !Zero.Equals(f(0, 1, Zero)))
{
var newRowPointers = RowPointers;
var newColumnIndices = new List<int>(ColumnIndices.Length);
var newValues = new List<T>(Values.Length);
int k = 0;
for (int row = 0; row < RowCount; row++)
{
newRowPointers[row] = newValues.Count;
for (int col = 0; col < ColumnCount; col++)
{
var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(row, col, Values[k++]) : f(row, col, Zero);
if (!Zero.Equals(item))
{
newValues.Add(item);
newColumnIndices.Add(col);
}
}
}
ColumnIndices = newColumnIndices.ToArray();
Values = newValues.ToArray();
newRowPointers[RowCount] = newValues.Count;
}
else
{
// we can safely do this in-place:
int nonZero = 0;
for (int row = 0; row < RowCount; row++)
{
var startIndex = RowPointers[row];
var endIndex = RowPointers[row + 1];
RowPointers[row] = nonZero;
for (var j = startIndex; j < endIndex; j++)
{
var item = f(row, ColumnIndices[j], Values[j]);
if (!Zero.Equals(item))
{
Values[nonZero] = item;
ColumnIndices[nonZero] = ColumnIndices[j];
nonZero++;
}
}
}
Array.Resize(ref ColumnIndices, nonZero);
Array.Resize(ref Values, nonZero);
RowPointers[RowCount] = nonZero;
}
}
internal override void MapToUnchecked<TU>(MatrixStorage<TU> target, Func<T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
{
var processZeros = forceMapZeros || !Zero.Equals(f(Zero));
@ -1212,62 +1318,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
{
if (forceMapZeros || !Zero.Equals(f(Zero)))
{
var newRowPointers = RowPointers;
var newColumnIndices = new List<int>(ColumnIndices.Length);
var newValues = new List<T>(Values.Length);
int k = 0;
for (int row = 0; row < RowCount; row++)
{
newRowPointers[row] = newValues.Count;
for (int col = 0; col < ColumnCount; col++)
{
var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(Values[k++]) : f(Zero);
if (!Zero.Equals(item))
{
newValues.Add(item);
newColumnIndices.Add(col);
}
}
}
ColumnIndices = newColumnIndices.ToArray();
Values = newValues.ToArray();
newRowPointers[RowCount] = newValues.Count;
}
else
{
// we can safely do this in-place:
int nonZero = 0;
for (int row = 0; row < RowCount; row++)
{
var startIndex = RowPointers[row];
var endIndex = RowPointers[row + 1];
RowPointers[row] = nonZero;
for (var j = startIndex; j < endIndex; j++)
{
var item = f(Values[j]);
if (!Zero.Equals(item))
{
Values[nonZero] = item;
ColumnIndices[nonZero] = ColumnIndices[j];
nonZero++;
}
}
}
Array.Resize(ref ColumnIndices, nonZero);
Array.Resize(ref Values, nonZero);
RowPointers[RowCount] = nonZero;
}
}
internal override void MapIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f, bool forceMapZeros = false, bool skipClearing = false)
{
var processZeros = forceMapZeros || !Zero.Equals(f(0, 0, Zero));
var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero));
var sparseTarget = target as SparseCompressedRowMatrixStorage<TU>;
if (sparseTarget != null)
@ -1359,56 +1412,213 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
internal override void MapSubMatrixIndexedToUnchecked<TU>(MatrixStorage<TU> target, Func<int, int, T, TU> f,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount,
bool forceMapZeros = false, bool skipClearing = false)
{
if (forceMapZeros || !Zero.Equals(f(0, 0, Zero)))
var sparseTarget = target as SparseCompressedRowMatrixStorage<TU>;
if (sparseTarget != null)
{
var newRowPointers = RowPointers;
var newColumnIndices = new List<int>(ColumnIndices.Length);
var newValues = new List<T>(Values.Length);
MapSubMatrixIndexedToUnchecked(sparseTarget, f, sourceRowIndex, targetRowIndex, rowCount, sourceColumnIndex, targetColumnIndex, columnCount, forceMapZeros, skipClearing);
return;
}
int k = 0;
for (int row = 0; row < RowCount; row++)
// FALL BACK
var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero));
if (!skipClearing && !processZeros)
{
target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
}
if (processZeros)
{
for (int sr = sourceRowIndex, tr = targetRowIndex; sr < sourceRowIndex + rowCount; sr++, tr++)
{
newRowPointers[row] = newValues.Count;
for (int col = 0; col < ColumnCount; col++)
var index = RowPointers[sr];
var endIndex = RowPointers[sr + 1];
// move forward to our sub-range
for (; ColumnIndices[index] < sourceColumnIndex && index < endIndex; index++)
{
var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(row, col, Values[k++]) : f(row, col, Zero);
if (!Zero.Equals(item))
}
for (int sc = sourceColumnIndex, tc = targetColumnIndex; sc < sourceColumnIndex + columnCount; sc++, tc++)
{
if (index < endIndex && sc == ColumnIndices[index])
{
newValues.Add(item);
newColumnIndices.Add(col);
target.At(tr, tc, f(tr, tc, Values[index]));
index = Math.Min(index + 1, endIndex);
}
else
{
target.At(tr, tc, f(tr, tc, Zero));
}
}
}
}
else
{
int columnOffset = targetColumnIndex - sourceColumnIndex;
for (int sr = sourceRowIndex, tr = targetRowIndex; sr < sourceRowIndex + rowCount; sr++, tr++)
{
var startIndex = RowPointers[sr];
var endIndex = RowPointers[sr + 1];
for (int k = startIndex; k < endIndex; k++)
{
// check if the column index is in the range
if ((ColumnIndices[k] >= sourceColumnIndex) && (ColumnIndices[k] < sourceColumnIndex + columnCount))
{
int tc = ColumnIndices[k] + columnOffset;
target.At(tr, tc, f(tr, tc, Values[k]));
}
}
}
}
}
ColumnIndices = newColumnIndices.ToArray();
Values = newValues.ToArray();
newRowPointers[RowCount] = newValues.Count;
void MapSubMatrixIndexedToUnchecked<TU>(SparseCompressedRowMatrixStorage<TU> target, Func<int, int, T, TU> f,
int sourceRowIndex, int targetRowIndex, int rowCount,
int sourceColumnIndex, int targetColumnIndex, int columnCount,
bool forceMapZeros = false, bool skipClearing = false)
where TU : struct, IEquatable<TU>, IFormattable
{
var processZeros = forceMapZeros || !Zero.Equals(f(0, 1, Zero));
if (!skipClearing && !processZeros)
{
target.Clear(targetRowIndex, rowCount, targetColumnIndex, columnCount);
}
var rowOffset = targetRowIndex - sourceRowIndex;
var columnOffset = targetColumnIndex - sourceColumnIndex;
var zero = Matrix<TU>.Zero;
// special case for empty target - much faster
if (target.ValueCount == 0)
{
var values = new List<TU>(ValueCount);
var columnIndices = new List<int>(ValueCount);
var rowPointers = target.RowPointers;
if (processZeros)
{
for (int sr = sourceRowIndex; sr < sourceRowIndex + rowCount; sr++)
{
int tr = sr + rowOffset;
rowPointers[tr] = values.Count;
var index = RowPointers[sr];
var endIndex = RowPointers[sr + 1];
// move forward to our sub-range
for (; ColumnIndices[index] < sourceColumnIndex && index < endIndex; index++)
{
}
for (int sc = sourceColumnIndex, tc = targetColumnIndex; sc < sourceColumnIndex + columnCount; sc++, tc++)
{
if (index < endIndex && sc == ColumnIndices[index])
{
TU item = f(tr, tc, Values[index]);
if (!zero.Equals(item))
{
values.Add(item);
columnIndices.Add(tc);
}
index = Math.Min(index + 1, endIndex);
}
else
{
TU item = f(tr, tc, Zero);
if (!zero.Equals(item))
{
values.Add(item);
columnIndices.Add(tc);
}
}
}
}
}
else
{
for (int sr = sourceRowIndex; sr < sourceRowIndex + rowCount; sr++)
{
int tr = sr + rowOffset;
rowPointers[tr] = values.Count;
var startIndex = RowPointers[sr];
var endIndex = RowPointers[sr + 1];
for (int k = startIndex; k < endIndex; k++)
{
// check if the column index is in the range
if ((ColumnIndices[k] >= sourceColumnIndex) && (ColumnIndices[k] < sourceColumnIndex + columnCount))
{
int tc = ColumnIndices[k] + columnOffset;
TU item = f(tr, tc, Values[k]);
if (!zero.Equals(item))
{
values.Add(item);
columnIndices.Add(tc);
}
}
}
}
}
for (int i = targetRowIndex + rowCount; i < rowPointers.Length; i++)
{
rowPointers[i] = values.Count;
}
target.RowPointers[target.RowCount] = values.Count;
target.Values = values.ToArray();
target.ColumnIndices = columnIndices.ToArray();
return;
}
// TODO: proper general sparse case - the following is essentially a fall back, not leveraging the target data structure
if (processZeros)
{
for (int sr = sourceRowIndex, tr = targetRowIndex; sr < sourceRowIndex + rowCount; sr++, tr++)
{
var index = RowPointers[sr];
var endIndex = RowPointers[sr + 1];
// move forward to our sub-range
for (; ColumnIndices[index] < sourceColumnIndex && index < endIndex; index++)
{
}
for (int sc = sourceColumnIndex, tc = targetColumnIndex; sc < sourceColumnIndex + columnCount; sc++, tc++)
{
if (index < endIndex && sc == ColumnIndices[index])
{
target.At(tr, tc, f(tr, tc, Values[index]));
index = Math.Min(index + 1, endIndex);
}
else
{
target.At(tr, tc, f(tr, tc, Zero));
}
}
}
}
else
{
// we can safely do this in-place:
int nonZero = 0;
for (int row = 0; row < RowCount; row++)
for (int sr = sourceRowIndex, tr = targetRowIndex; sr < sourceRowIndex + rowCount; sr++, tr++)
{
var startIndex = RowPointers[row];
var endIndex = RowPointers[row + 1];
RowPointers[row] = nonZero;
for (var j = startIndex; j < endIndex; j++)
var startIndex = RowPointers[sr];
var endIndex = RowPointers[sr + 1];
for (int k = startIndex; k < endIndex; k++)
{
var item = f(row, ColumnIndices[j], Values[j]);
if (!Zero.Equals(item))
// check if the column index is in the range
if ((ColumnIndices[k] >= sourceColumnIndex) && (ColumnIndices[k] < sourceColumnIndex + columnCount))
{
Values[nonZero] = item;
ColumnIndices[nonZero] = ColumnIndices[j];
nonZero++;
int tc = ColumnIndices[k] + columnOffset;
target.At(tr, tc, f(tr, tc, Values[k]));
}
}
}
Array.Resize(ref ColumnIndices, nonZero);
Array.Resize(ref Values, nonZero);
RowPointers[RowCount] = nonZero;
}
}
}

260
src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Map.cs

@ -0,0 +1,260 @@
// <copyright file="MatrixStructureTheory.Map.cs" company="Math.NET">
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
// http://mathnetnumerics.codeplex.com
//
// Copyright (c) 2009-2014 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using System.Linq;
using MathNet.Numerics.LinearAlgebra;
using NUnit.Framework;
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests
{
partial class MatrixStructureTheory<T>
{
[Theory]
public void CanMap(Matrix<T> matrix)
{
Matrix<T> a = matrix.Map(x => x, false);
Assert.That(a, Is.EqualTo(matrix));
Assert.That(a.Storage.IsDense, Is.EqualTo(matrix.Storage.IsDense));
Assert.That(a.Storage.IsFullyMutable, Is.EqualTo(matrix.Storage.IsFullyMutable));
T one = Matrix<T>.Build.One;
Assert.That(matrix.Map(x => x, true), Is.EqualTo(matrix));
Assert.That(matrix.Map(x => one, true), Is.EqualTo(Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one)));
// Map into existing - we skip zeros, but existing values must still be reset to zero
var dense = Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one);
matrix.Map(x => x, dense, false);
Assert.That(dense, Is.EqualTo(matrix));
// Map into self, without using the proper MapInplace method:
var copy = matrix.Clone();
copy.Map(x => x, copy, false);
Assert.That(copy, Is.EqualTo(matrix));
if (matrix.Storage.IsFullyMutable)
{
copy.Map(x => one, copy, false);
Assert.That(copy, Is.EqualTo(Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one)));
}
}
[Theory]
public void CanMapIndexed(Matrix<T> matrix)
{
Matrix<T> a = matrix.MapIndexed((i, j, x) =>
{
if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x));
return x;
}, false);
Assert.That(a, Is.EqualTo(matrix));
Assert.That(a.Storage.IsDense, Is.EqualTo(matrix.Storage.IsDense));
Assert.That(a.Storage.IsFullyMutable, Is.EqualTo(matrix.Storage.IsFullyMutable));
T one = Matrix<T>.Build.One;
var d = matrix.MapIndexed((i, j, x) => i == j ? one : x, false);
Assert.That(d.Diagonal().All(x => one.Equals(x) || Zero.Equals(x)), Is.True);
Assert.That(d.EnumerateIndexed().All(z => (z.Item1 == z.Item2) || (matrix.At(z.Item1, z.Item2).Equals(z.Item3))));
Assert.That(matrix.MapIndexed((i, j, x) => x, true), Is.EqualTo(matrix));
Assert.That(matrix.MapIndexed((i, j, x) => one, true), Is.EqualTo(Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one)));
// Map into existing - we skip zeros, but existing values must still be reset to zero
var dense = Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one);
matrix.MapIndexed((i, j, x) => x, dense, false);
Assert.That(dense, Is.EqualTo(matrix));
// Map into self, without using the proper MapInplace method:
var copy = matrix.Clone();
copy.MapIndexed((i, j, x) => x, copy, false);
Assert.That(copy, Is.EqualTo(matrix));
if (matrix.Storage.IsFullyMutable)
{
copy.MapIndexed((i, j, x) => one, copy, false);
Assert.That(copy, Is.EqualTo(Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one)));
}
}
[Theory]
public void CanMapInplace(Matrix<T> matrix)
{
var a = matrix.Clone();
a.MapInplace(x => x, false);
Assert.That(a, Is.EqualTo(matrix));
if (matrix.Storage.IsFullyMutable)
{
a.MapInplace(x => x, true);
Assert.That(a, Is.EqualTo(matrix));
T one = Matrix<T>.Build.One;
a.MapInplace(x => one, true);
Assert.That(a, Is.EqualTo(Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one)));
}
}
[Theory]
public void CanMapIndexedInplace(Matrix<T> matrix)
{
var a = matrix.Clone();
a.MapIndexedInplace((i, j, x) =>
{
if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x));
return x;
}, false);
Assert.That(a, Is.EqualTo(matrix));
if (matrix.Storage.IsFullyMutable)
{
a.MapIndexedInplace((i, j, x) => x, true);
Assert.That(a, Is.EqualTo(matrix));
T one = Matrix<T>.Build.One;
a.MapIndexedInplace((i, j, x) => i == j ? one : x, false);
Assert.That(a.Diagonal().All(x => one.Equals(x) || Zero.Equals(x)), Is.True);
Assert.That(a.EnumerateIndexed().All(z => (z.Item1 == z.Item2) || (matrix.At(z.Item1, z.Item2).Equals(z.Item3))));
a.MapIndexedInplace((i, j, x) => one, true);
Assert.That(a, Is.EqualTo(Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one)));
}
}
[Theory]
public void CanMapSubMatrixToSame(Matrix<T> matrix)
{
T one = Matrix<T>.Build.One;
// Full Range - not forced
Matrix<T> target = Matrix<T>.Build.SameAs(matrix);
matrix.Storage.MapSubMatrixIndexedTo(target.Storage, (i, j, x) => x, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false);
Assert.That(target, Is.EqualTo(matrix), "Full Range - not forced");
matrix.Storage.MapSubMatrixIndexedTo(target.Storage, (i, j, x) => Zero.Equals(x) ? Zero : one, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false);
Assert.That(target.Enumerate().All(x => Zero.Equals(x) || one.Equals(x)), Is.True);
if (matrix.Storage.IsFullyMutable)
{
// Full Range - forced
target = Matrix<T>.Build.SameAs(matrix);
matrix.Storage.MapSubMatrixIndexedTo(target.Storage, (i, j, x) => x, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, true, false);
Assert.That(target, Is.EqualTo(matrix), "Full Range - forced");
matrix.Storage.MapSubMatrixIndexedTo(target.Storage, (i, j, x) => Zero.Equals(x) ? Zero : one, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, true, false);
Assert.That(target.Enumerate().All(x => Zero.Equals(x) || one.Equals(x)), Is.True);
}
}
[Theory]
public void CanMapSubMatrixToDense(Matrix<T> matrix)
{
T one = Matrix<T>.Build.One;
// Full Range - not forced
Matrix<T> dense = Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one);
matrix.Storage.MapSubMatrixIndexedTo(dense.Storage, (i, j, x) =>
{
if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x));
return x;
}, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false);
Assert.That(dense, Is.EqualTo(matrix), "Full Range - not forced");
// Full Range - forced
dense = Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount, one);
matrix.Storage.MapSubMatrixIndexedTo(dense.Storage, (i, j, x) =>
{
if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x));
return x;
}, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, true, false);
Assert.That(dense, Is.EqualTo(matrix), "Full Range - forced");
// Sub Range - not forced - all except first column padded into 1-border
dense = Matrix<T>.Build.Dense(matrix.RowCount + 2, matrix.ColumnCount + 1, one);
matrix.Storage.MapSubMatrixIndexedTo(dense.Storage, (i, j, x) => x, 0, 1, matrix.RowCount, 1, 1, matrix.ColumnCount - 1, false, false);
Assert.That(dense.SubMatrix(1, dense.RowCount - 2, 1, dense.ColumnCount - 2),
Is.EqualTo(matrix.SubMatrix(0, matrix.RowCount, 1, matrix.ColumnCount - 1)), "Sub Range - not forced - range");
dense.SetSubMatrix(1, 0, matrix.RowCount, 1, 0, matrix.ColumnCount - 1, Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount - 1, one));
Assert.That(dense.Enumerate().All(one.Equals), Is.True);
// Sub Range - forced - all except first row padded into 1-border
dense = Matrix<T>.Build.Dense(matrix.RowCount + 1, matrix.ColumnCount + 2, one);
matrix.Storage.MapSubMatrixIndexedTo(dense.Storage, (i, j, x) => x, 1, 1, matrix.RowCount - 1, 0, 1, matrix.ColumnCount, true, false);
Assert.That(dense.SubMatrix(1, dense.RowCount - 2, 1, dense.ColumnCount - 2),
Is.EqualTo(matrix.SubMatrix(1, matrix.RowCount - 1, 0, matrix.ColumnCount)), "Sub Range - forced - range");
dense.SetSubMatrix(1, 0, matrix.RowCount - 1, 1, 0, matrix.ColumnCount, Matrix<T>.Build.Dense(matrix.RowCount - 1, matrix.ColumnCount, one));
Assert.That(dense.Enumerate().All(one.Equals), Is.True);
}
[Theory]
public void CanMapSubMatrixToSparse(Matrix<T> matrix)
{
T one = Matrix<T>.Build.One;
// Full Range - filled, not forced
Matrix<T> sparse = Matrix<T>.Build.Sparse(matrix.RowCount, matrix.ColumnCount, one);
matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) =>
{
if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x));
return x;
}, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false);
Assert.That(sparse, Is.EqualTo(matrix), "Full Range - filled, not forced");
// Full Range - empty, not forced
sparse = Matrix<T>.Build.Sparse(matrix.RowCount, matrix.ColumnCount);
matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) =>
{
if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x));
return x;
}, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, false, false);
Assert.That(sparse, Is.EqualTo(matrix), "Full Range - empty, not forced");
// Full Range - filled, forced
sparse = Matrix<T>.Build.Sparse(matrix.RowCount, matrix.ColumnCount, one);
matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) =>
{
if (i != 0 || j != 1) Assert.That(matrix.At(i, j), Is.EqualTo(x));
return x;
}, 0, 0, matrix.RowCount, 0, 0, matrix.ColumnCount, true, false);
Assert.That(sparse, Is.EqualTo(matrix), "Full Range - filled, forced");
// Sub Range - filled, not forced - all except first column padded into 1-border
sparse = Matrix<T>.Build.Sparse(matrix.RowCount + 2, matrix.ColumnCount + 1, one);
matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) => x, 0, 1, matrix.RowCount, 1, 1, matrix.ColumnCount - 1, false, false);
Assert.That(sparse.SubMatrix(1, sparse.RowCount - 2, 1, sparse.ColumnCount - 2),
Is.EqualTo(matrix.SubMatrix(0, matrix.RowCount, 1, matrix.ColumnCount - 1)), "Sub Range - filled, not forced - range");
sparse.SetSubMatrix(1, 0, matrix.RowCount, 1, 0, matrix.ColumnCount - 1, Matrix<T>.Build.Dense(matrix.RowCount, matrix.ColumnCount - 1, one));
Assert.That(sparse.Enumerate().All(one.Equals), Is.True);
// Sub Range - filled, forced - all except first row padded into 1-border
sparse = Matrix<T>.Build.Sparse(matrix.RowCount + 1, matrix.ColumnCount + 2, one);
matrix.Storage.MapSubMatrixIndexedTo(sparse.Storage, (i, j, x) => x, 1, 1, matrix.RowCount - 1, 0, 1, matrix.ColumnCount, true, false);
Assert.That(sparse.SubMatrix(1, sparse.RowCount - 2, 1, sparse.ColumnCount - 2),
Is.EqualTo(matrix.SubMatrix(1, matrix.RowCount - 1, 0, matrix.ColumnCount)), "Sub Range - filled, forced - range");
sparse.SetSubMatrix(1, 0, matrix.RowCount - 1, 1, 0, matrix.ColumnCount, Matrix<T>.Build.Dense(matrix.RowCount - 1, matrix.ColumnCount, one));
Assert.That(sparse.Enumerate().All(one.Equals), Is.True);
}
}
}

4
src/UnitTests/LinearAlgebraTests/MatrixStructureTheory.Reform.cs

@ -28,10 +28,10 @@
// OTHER DEALINGS IN THE SOFTWARE.
// </copyright>
using MathNet.Numerics.LinearAlgebra;
using NUnit.Framework;
using System;
using System.Linq;
using MathNet.Numerics.LinearAlgebra;
using NUnit.Framework;
namespace MathNet.Numerics.UnitTests.LinearAlgebraTests
{

1
src/UnitTests/UnitTests.csproj

@ -305,6 +305,7 @@
<Compile Include="LinearAlgebraTests\Double\VectorTests.Norm.cs" />
<Compile Include="LinearAlgebraTests\MatrixStructureTheory.Access.cs" />
<Compile Include="LinearAlgebraTests\MatrixStructureTheory.cs" />
<Compile Include="LinearAlgebraTests\MatrixStructureTheory.Map.cs" />
<Compile Include="LinearAlgebraTests\MatrixStructureTheory.Reform.cs" />
<Compile Include="LinearAlgebraTests\Single\DenseMatrixTests.cs" />
<Compile Include="LinearAlgebraTests\Single\DenseVectorArithmeticTheory.cs" />

Loading…
Cancel
Save