Browse Source

LA: Migrate sparse CSR format to more common row pointer convention

pull/163/head
Christoph Ruegg 13 years ago
parent
commit
afec8f3c2f
  1. 142
      src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
  2. 142
      src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
  3. 11
      src/Numerics/LinearAlgebra/Double/Solvers/MILU0Preconditioner.cs
  4. 141
      src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
  5. 133
      src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
  6. 173
      src/Numerics/LinearAlgebra/Storage/SparseCompressedRowMatrixStorage.cs

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

@ -379,13 +379,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row >= columnIndices[j])
{
@ -446,13 +444,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row <= columnIndices[j])
{
@ -514,13 +510,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row > columnIndices[j])
{
@ -582,13 +576,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row < columnIndices[j])
{
@ -607,25 +599,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var ret = new SparseCompressedRowMatrixStorage<Complex>(ColumnCount, RowCount)
{
ColumnIndices = new int[valueCount],
Values = new Complex[valueCount]
ColumnIndices = new int[_storage.ValueCount],
Values = new Complex[_storage.ValueCount]
};
// Do an 'inverse' CopyTo iterate over the rows
for (var i = 0; i < rowPointers.Length; i++)
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] : NonZerosCount;
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;
}
@ -645,16 +633,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var aat = (SparseCompressedRowMatrixStorage<Complex>) (this*ConjugateTranspose()).Storage;
var norm = 0d;
for (var i = 0; i < aat.RowPointers.Length; i++)
for (var i = 0; i < aat.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = aat.RowPointers[i];
var endIndex = i < aat.RowPointers.Length - 1 ? aat.RowPointers[i + 1] : aat.ValueCount;
var endIndex = aat.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;
}
@ -677,19 +662,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
{
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var norm = 0d;
for (var i = 0; i < rowPointers.Length; i++)
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;
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;
}
@ -716,21 +697,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </exception>
public static SparseMatrix Identity(int order)
{
var m = new SparseCompressedRowMatrixStorage<Complex>(order, order)
var storage = new SparseCompressedRowMatrixStorage<Complex>(order, order)
{
ValueCount = order,
Values = new Complex[order],
ColumnIndices = new int[order]
};
for (var i = 0; i < order; i++)
{
m.Values[i] = 1d;
m.ColumnIndices[i] = i;
m.RowPointers[i] = i;
storage.Values[i] = 1d;
storage.ColumnIndices[i] = i;
storage.RowPointers[i] = i;
}
return new SparseMatrix(m);
storage.RowPointers[order] = order;
return new SparseMatrix(storage);
}
#endregion
@ -781,11 +763,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = leftStorage.RowPointers[i];
var endIndex = i < leftStorage.RowPointers.Length - 1 ? leftStorage.RowPointers[i + 1] : leftStorage.ValueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = leftStorage.RowPointers[i + 1];
for (var j = leftStorage.RowPointers[i]; j < endIndex; j++)
{
var columnIndex = leftStorage.ColumnIndices[j];
var resVal = leftStorage.Values[j] + result.At(i, columnIndex);
@ -824,11 +803,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
{
for (var i = 0; i < otherStorage.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = otherStorage.RowPointers[i];
var endIndex = i < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[i + 1] : otherStorage.ValueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = otherStorage.RowPointers[i + 1];
for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
{
var columnIndex = otherStorage.ColumnIndices[j];
var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
@ -848,15 +824,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
var columnIndex = columnIndices[j];
var resVal = sparseResult.At(i, columnIndex) + values[j];
@ -935,13 +907,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < RowCount; row++)
{
// Get the begin / end index for the current row
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
var endIndex = rowPointers[row + 1];
if (startIndex == endIndex)
{
continue;
@ -973,13 +944,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < RowCount; row++)
{
// Get the begin / end index for the current row
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
var endIndex = rowPointers[row + 1];
if (startIndex == endIndex)
{
continue;
@ -1015,15 +985,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var otherStorage = otherSparse._storage;
for (var j = 0; j < RowCount; j++)
{
// Get the begin / end index for the row
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = j < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[j + 1] : otherStorage.ValueCount;
var endIndexOther = otherStorage.RowPointers[j + 1];
if (startIndexOther == endIndexOther)
{
continue;
@ -1032,9 +1001,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
for (var i = 0; i < RowCount; i++)
{
// Multiply row of matrix A on row of matrix B
// Get the begin / end index for the row
var startIndexThis = rowPointers[i];
var endIndexThis = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount;
var endIndexThis = rowPointers[i + 1];
if (startIndexThis == endIndexThis)
{
continue;
@ -1077,15 +1047,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
var resVal = values[j]*other.At(i, columnIndices[j]);
if (!resVal.IsZero())
@ -1108,15 +1074,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
if (!values[j].IsZero())
{
@ -1146,15 +1108,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
if (!values[j].IsZero())
{
@ -1178,7 +1136,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
// todo: we might be able to speed this up by caching one half of the matrix
var rowPointers = _storage.RowPointers;
for (var row = 0; row < RowCount - 1; row++)
for (var row = 0; row < RowCount; row++)
{
var start = rowPointers[row];
var end = rowPointers[row + 1];
@ -1194,16 +1152,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
}
}
var lastRow = rowPointers.Length - 1;
if (rowPointers[lastRow] < NonZerosCount)
{
if (!CheckIfOppositesAreEqual(rowPointers[lastRow], _storage.ValueCount, lastRow))
{
return false;
}
}
return true;
}
}

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

@ -374,13 +374,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row >= columnIndices[j])
{
@ -441,13 +439,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row <= columnIndices[j])
{
@ -509,13 +505,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row > columnIndices[j])
{
@ -577,13 +571,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row < columnIndices[j])
{
@ -602,25 +594,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var ret = new SparseCompressedRowMatrixStorage<Complex32>(ColumnCount, RowCount)
{
ColumnIndices = new int[valueCount],
Values = new Complex32[valueCount]
ColumnIndices = new int[_storage.ValueCount],
Values = new Complex32[_storage.ValueCount]
};
// Do an 'inverse' CopyTo iterate over the rows
for (var i = 0; i < rowPointers.Length; i++)
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] : NonZerosCount;
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;
}
@ -640,16 +628,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var aat = (SparseCompressedRowMatrixStorage<Complex32>) (this*ConjugateTranspose()).Storage;
var norm = 0f;
for (var i = 0; i < aat.RowPointers.Length; i++)
for (var i = 0; i < aat.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = aat.RowPointers[i];
var endIndex = i < aat.RowPointers.Length - 1 ? aat.RowPointers[i + 1] : aat.ValueCount;
var endIndex = aat.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;
}
@ -672,19 +657,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
{
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var norm = 0f;
for (var i = 0; i < rowPointers.Length; i++)
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;
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;
}
@ -711,21 +692,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </exception>
public static SparseMatrix Identity(int order)
{
var m = new SparseCompressedRowMatrixStorage<Complex32>(order, order)
var storage = new SparseCompressedRowMatrixStorage<Complex32>(order, order)
{
ValueCount = order,
Values = new Complex32[order],
ColumnIndices = new int[order]
};
for (var i = 0; i < order; i++)
{
m.Values[i] = 1f;
m.ColumnIndices[i] = i;
m.RowPointers[i] = i;
storage.Values[i] = 1f;
storage.ColumnIndices[i] = i;
storage.RowPointers[i] = i;
}
return new SparseMatrix(m);
storage.RowPointers[order] = order;
return new SparseMatrix(storage);
}
#endregion
@ -776,11 +758,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = leftStorage.RowPointers[i];
var endIndex = i < leftStorage.RowPointers.Length - 1 ? leftStorage.RowPointers[i + 1] : leftStorage.ValueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = leftStorage.RowPointers[i + 1];
for (var j = leftStorage.RowPointers[i]; j < endIndex; j++)
{
var columnIndex = leftStorage.ColumnIndices[j];
var resVal = leftStorage.Values[j] + result.At(i, columnIndex);
@ -818,11 +797,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
{
for (var i = 0; i < otherStorage.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = otherStorage.RowPointers[i];
var endIndex = i < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[i + 1] : otherStorage.ValueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = otherStorage.RowPointers[i + 1];
for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
{
var columnIndex = otherStorage.ColumnIndices[j];
var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
@ -842,15 +818,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
var columnIndex = columnIndices[j];
var resVal = sparseResult.At(i, columnIndex) + values[j];
@ -929,13 +901,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < RowCount; row++)
{
// Get the begin / end index for the current row
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
var endIndex = rowPointers[row + 1];
if (startIndex == endIndex)
{
continue;
@ -967,13 +938,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < RowCount; row++)
{
// Get the begin / end index for the current row
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
var endIndex = rowPointers[row + 1];
if (startIndex == endIndex)
{
continue;
@ -1009,15 +979,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var otherStorage = otherSparse._storage;
for (var j = 0; j < RowCount; j++)
{
// Get the begin / end index for the row
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = j < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[j + 1] : otherStorage.ValueCount;
var endIndexOther = otherStorage.RowPointers[j + 1];
if (startIndexOther == endIndexOther)
{
continue;
@ -1026,9 +995,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
for (var i = 0; i < RowCount; i++)
{
// Multiply row of matrix A on row of matrix B
// Get the begin / end index for the row
var startIndexThis = rowPointers[i];
var endIndexThis = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount;
var endIndexThis = rowPointers[i + 1];
if (startIndexThis == endIndexThis)
{
continue;
@ -1071,15 +1041,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
var resVal = values[j]*other.At(i, columnIndices[j]);
if (!resVal.IsZero())
@ -1102,15 +1068,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
if (!values[j].IsZero())
{
@ -1140,15 +1102,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
if (!values[j].IsZero())
{
@ -1172,7 +1130,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
// todo: we might be able to speed this up by caching one half of the matrix
var rowPointers = _storage.RowPointers;
for (var row = 0; row < RowCount - 1; row++)
for (var row = 0; row < RowCount; row++)
{
var start = rowPointers[row];
var end = rowPointers[row + 1];
@ -1188,16 +1146,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
}
}
var lastRow = rowPointers.Length - 1;
if (rowPointers[lastRow] < NonZerosCount)
{
if (!CheckIfOppositesAreEqual(rowPointers[lastRow], _storage.ValueCount, lastRow))
{
return false;
}
}
return true;
}
}

11
src/Numerics/LinearAlgebra/Double/Solvers/MILU0Preconditioner.cs

@ -116,15 +116,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers
// Original matrix compressed sparse row storage.
int[] ja = csr.ColumnIndices;
double[] a = csr.Values;
int[] ia = csr.RowPointers;
// Make row pointers conform to standard CSR definition (ie. length n+1)
// TODO: We plan to migrate to this definition everywhere
var ia = new int[n + 1];
Array.Copy(csr.RowPointers, ia, n);
ia[n] = csr.ValueCount;
_alu = new double[csr.ValueCount + 1];
_jlu = new int[csr.ValueCount + 1];
_alu = new double[ia.Length];
_jlu = new int[ia.Length];
_diag = new int[n];
int code = Compute(n, a, ja, ia, _alu, _jlu, _diag, UseModified);

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

@ -372,13 +372,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row >= columnIndices[j])
{
@ -439,13 +437,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row <= columnIndices[j])
{
@ -507,13 +503,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row > columnIndices[j])
{
@ -575,13 +569,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row < columnIndices[j])
{
@ -600,22 +592,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var ret = new SparseCompressedRowMatrixStorage<double>(ColumnCount, RowCount)
{
ColumnIndices = new int[valueCount],
Values = new double[valueCount]
ColumnIndices = new int[_storage.ValueCount],
Values = new double[_storage.ValueCount]
};
// Do an 'inverse' CopyTo iterate over the rows
for (var i = 0; i < rowPointers.Length; i++)
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;
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
@ -638,13 +627,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var aat = (SparseCompressedRowMatrixStorage<double>) (this*Transpose()).Storage;
var norm = 0d;
for (var i = 0; i < aat.RowPointers.Length; i++)
for (var i = 0; i < aat.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = aat.RowPointers[i];
var endIndex = i < aat.RowPointers.Length - 1 ? aat.RowPointers[i + 1] : aat.ValueCount;
var endIndex = aat.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
@ -669,16 +656,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var norm = 0d;
for (var i = 0; i < rowPointers.Length; i++)
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;
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
@ -708,21 +692,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </exception>
public static SparseMatrix Identity(int order)
{
var m = new SparseCompressedRowMatrixStorage<double>(order, order)
var storage = new SparseCompressedRowMatrixStorage<double>(order, order)
{
ValueCount = order,
Values = new double[order],
ColumnIndices = new int[order]
};
for (var i = 0; i < order; i++)
{
m.Values[i] = 1d;
m.ColumnIndices[i] = i;
m.RowPointers[i] = i;
storage.Values[i] = 1d;
storage.ColumnIndices[i] = i;
storage.RowPointers[i] = i;
}
return new SparseMatrix(m);
storage.RowPointers[order] = order;
return new SparseMatrix(storage);
}
#endregion
@ -773,11 +758,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = leftStorage.RowPointers[i];
var endIndex = i < leftStorage.RowPointers.Length - 1 ? leftStorage.RowPointers[i + 1] : leftStorage.ValueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = leftStorage.RowPointers[i + 1];
for (var j = leftStorage.RowPointers[i]; j < endIndex; j++)
{
var columnIndex = leftStorage.ColumnIndices[j];
var resVal = leftStorage.Values[j] + result.At(i, columnIndex);
@ -816,11 +798,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double
{
for (var i = 0; i < otherStorage.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = otherStorage.RowPointers[i];
var endIndex = i < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[i + 1] : otherStorage.ValueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = otherStorage.RowPointers[i + 1];
for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
{
var columnIndex = otherStorage.ColumnIndices[j];
var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
@ -840,15 +819,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
var columnIndex = columnIndices[j];
var resVal = sparseResult.At(i, columnIndex) + values[j];
@ -871,7 +846,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return;
}
if (scalar == 0.0 || _storage.ValueCount == 0)
if (scalar == 0.0 || NonZerosCount == 0)
{
result.Clear();
return;
@ -927,13 +902,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < RowCount; row++)
{
// Get the begin / end index for the current row
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
var endIndex = rowPointers[row + 1];
if (startIndex == endIndex)
{
continue;
@ -965,13 +939,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < RowCount; row++)
{
// Get the begin / end index for the current row
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
var endIndex = rowPointers[row + 1];
if (startIndex == endIndex)
{
continue;
@ -1007,15 +980,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var otherStorage = otherSparse._storage;
for (var j = 0; j < RowCount; j++)
{
// Get the begin / end index for the row
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = j < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[j + 1] : otherStorage.ValueCount;
var endIndexOther = otherStorage.RowPointers[j + 1];
if (startIndexOther == endIndexOther)
{
continue;
@ -1023,10 +995,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
for (var i = 0; i < RowCount; i++)
{
// Multiply row of matrix A on row of matrix B
// Get the begin / end index for the row
var startIndexThis = rowPointers[i];
var endIndexThis = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount;
var endIndexThis = rowPointers[i + 1];
if (startIndexThis == endIndexThis)
{
continue;
@ -1069,15 +1040,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
var resVal = values[j]*other.At(i, columnIndices[j]);
if (resVal != 0d)
@ -1100,15 +1067,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
if (values[j] != 0d)
{
@ -1138,15 +1101,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
if (values[j] != 0d)
{
@ -1196,7 +1155,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
// todo: we might be able to speed this up by caching one half of the matrix
var rowPointers = _storage.RowPointers;
for (var row = 0; row < RowCount - 1; row++)
for (var row = 0; row < RowCount; row++)
{
var start = rowPointers[row];
var end = rowPointers[row + 1];
@ -1212,16 +1171,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
}
}
var lastRow = rowPointers.Length - 1;
if (rowPointers[lastRow] < _storage.ValueCount)
{
if (!CheckIfOppositesAreEqual(rowPointers[lastRow], _storage.ValueCount, lastRow))
{
return false;
}
}
return true;
}
}

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

@ -372,13 +372,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row >= columnIndices[j])
{
@ -439,13 +437,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row <= columnIndices[j])
{
@ -507,13 +503,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row > columnIndices[j])
{
@ -575,13 +569,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < result.RowCount; row++)
{
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = rowPointers[row + 1];
for (var j = rowPointers[row]; j < endIndex; j++)
{
if (row < columnIndices[j])
{
@ -600,20 +592,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var ret = new SparseCompressedRowMatrixStorage<float>(ColumnCount, RowCount)
{
ColumnIndices = new int[valueCount],
Values = new float[valueCount]
ColumnIndices = new int[_storage.ValueCount],
Values = new float[_storage.ValueCount]
};
// Do an 'inverse' CopyTo iterate over the rows
for (var i = 0; i < rowPointers.Length; i++)
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;
var endIndex = rowPointers[i + 1];
// Get the values for the current row
if (startIndex == endIndex)
@ -638,11 +629,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var aat = (SparseCompressedRowMatrixStorage<float>) (this*Transpose()).Storage;
var norm = 0f;
for (var i = 0; i < aat.RowPointers.Length; i++)
for (var i = 0; i < aat.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = aat.RowPointers[i];
var endIndex = i < aat.RowPointers.Length - 1 ? aat.RowPointers[i + 1] : aat.ValueCount;
var endIndex = aat.RowPointers[i + 1];
// Get the values for the current row
if (startIndex == endIndex)
@ -669,14 +660,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
{
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var norm = 0f;
for (var i = 0; i < rowPointers.Length; i++)
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;
var endIndex = rowPointers[i + 1];
// Get the values for the current row
if (startIndex == endIndex)
@ -708,21 +698,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </exception>
public static SparseMatrix Identity(int order)
{
var mStorage = new SparseCompressedRowMatrixStorage<float>(order, order)
var storage = new SparseCompressedRowMatrixStorage<float>(order, order)
{
ValueCount = order,
Values = new float[order],
ColumnIndices = new int[order]
};
for (var i = 0; i < order; i++)
{
mStorage.Values[i] = 1f;
mStorage.ColumnIndices[i] = i;
mStorage.RowPointers[i] = i;
storage.Values[i] = 1f;
storage.ColumnIndices[i] = i;
storage.RowPointers[i] = i;
}
return new SparseMatrix(mStorage);
storage.RowPointers[order] = order;
return new SparseMatrix(storage);
}
#endregion
@ -773,11 +764,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = leftStorage.RowPointers[i];
var endIndex = i < leftStorage.RowPointers.Length - 1 ? leftStorage.RowPointers[i + 1] : leftStorage.ValueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = leftStorage.RowPointers[i + 1];
for (var j = leftStorage.RowPointers[i]; j < endIndex; j++)
{
var columnIndex = leftStorage.ColumnIndices[j];
var resVal = leftStorage.Values[j] + result.At(i, columnIndex);
@ -815,11 +803,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single
{
for (var i = 0; i < otherStorage.RowCount; i++)
{
// Get the begin / end index for the current row
var startIndex = otherStorage.RowPointers[i];
var endIndex = i < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[i + 1] : otherStorage.ValueCount;
for (var j = startIndex; j < endIndex; j++)
var endIndex = otherStorage.RowPointers[i + 1];
for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
{
var columnIndex = otherStorage.ColumnIndices[j];
var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
@ -839,15 +824,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
var columnIndex = columnIndices[j];
var resVal = sparseResult.At(i, columnIndex) + values[j];
@ -926,13 +907,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < RowCount; row++)
{
// Get the begin / end index for the current row
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
var endIndex = rowPointers[row + 1];
if (startIndex == endIndex)
{
continue;
@ -964,13 +944,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
for (var row = 0; row < RowCount; row++)
{
// Get the begin / end index for the current row
var startIndex = rowPointers[row];
var endIndex = row < rowPointers.Length - 1 ? rowPointers[row + 1] : valueCount;
var endIndex = rowPointers[row + 1];
if (startIndex == endIndex)
{
continue;
@ -1006,15 +985,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var valueCount = _storage.ValueCount;
var otherStorage = otherSparse._storage;
for (var j = 0; j < RowCount; j++)
{
// Get the begin / end index for the row
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = j < otherStorage.RowPointers.Length - 1 ? otherStorage.RowPointers[j + 1] : otherStorage.ValueCount;
var endIndexOther = otherStorage.RowPointers[j + 1];
if (startIndexOther == endIndexOther)
{
continue;
@ -1023,9 +1001,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
for (var i = 0; i < RowCount; i++)
{
// Multiply row of matrix A on row of matrix B
// Get the begin / end index for the row
var startIndexThis = rowPointers[i];
var endIndexThis = i < rowPointers.Length - 1 ? rowPointers[i + 1] : valueCount;
var endIndexThis = rowPointers[i + 1];
if (startIndexThis == endIndexThis)
{
continue;
@ -1068,15 +1047,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
var resVal = values[j]*other.At(i, columnIndices[j]);
if (resVal != 0f)
@ -1099,15 +1074,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
if (values[j] != 0f)
{
@ -1137,15 +1108,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single
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++)
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
{
if (values[j] != 0f)
{
@ -1195,7 +1162,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
// todo: we might be able to speed this up by caching one half of the matrix
var rowPointers = _storage.RowPointers;
for (var row = 0; row < RowCount - 1; row++)
for (var row = 0; row < RowCount; row++)
{
var start = rowPointers[row];
var end = rowPointers[row + 1];
@ -1211,16 +1178,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
}
}
var lastRow = rowPointers.Length - 1;
if (rowPointers[lastRow] < NonZerosCount)
{
if (!CheckIfOppositesAreEqual(rowPointers[lastRow], NonZerosCount, lastRow))
{
return false;
}
}
return true;
}
}

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

@ -42,14 +42,16 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// [ruegg] public fields are OK here
/// <summary>
/// The array containing the row indices of the existing rows. Element "j" of the array gives the index of the
/// element in the <see cref="Values"/> array that is first non-zero element in a row "j"
/// The array containing the row indices of the existing rows. Element "i" of the array gives the index of the
/// element in the <see cref="Values"/> array that is first non-zero element in a row "i".
/// The last value is equal to ValueCount, so that the number of non-zero entries in row "i" is always
/// given by RowPointers[i+i] - RowPointers[i]. This array thus has length RowCount+1.
/// </summary>
public readonly int[] RowPointers;
/// <summary>
/// An array containing the column indices of the non-zero values. Element "I" of the array
/// is the number of the column in matrix that contains the I-th value in the <see cref="Values"/> array.
/// An array containing the column indices of the non-zero values. Element "j" of the array
/// is the number of the column in matrix that contains the j-th value in the <see cref="Values"/> array.
/// </summary>
public int[] ColumnIndices;
@ -63,15 +65,17 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
/// Gets the number of non zero elements in the matrix.
/// </summary>
/// <value>The number of non zero elements.</value>
public int ValueCount;
public int ValueCount
{
get { return RowPointers[RowCount]; }
}
internal SparseCompressedRowMatrixStorage(int rows, int columns)
: base(rows, columns)
{
RowPointers = new int[rows];
RowPointers = new int[rows + 1];
ColumnIndices = new int[0];
Values = new T[0];
ValueCount = 0;
}
/// <summary>
@ -152,34 +156,32 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
index = ~index;
var valueCount = RowPointers[RowPointers.Length - 1];
// Check if the storage needs to be increased
if ((ValueCount == Values.Length) && (ValueCount < ((long)RowCount * ColumnCount)))
if ((valueCount == Values.Length) && (valueCount < ((long)RowCount * ColumnCount)))
{
// Value array is completely full so we increase the size
// Determine the increase in size. We will not grow beyond the size of the matrix
var size = Math.Min(Values.Length + GrowthSize(), (long)RowCount * ColumnCount);
var size = Math.Min(Values.Length + GrowthSize(), (long) RowCount*ColumnCount);
if (size > int.MaxValue)
{
throw new NotSupportedException(Resources.TooManyElements);
}
Array.Resize(ref Values, (int)size);
Array.Resize(ref ColumnIndices, (int)size);
Array.Resize(ref Values, (int) size);
Array.Resize(ref ColumnIndices, (int) size);
}
// Move all values (with a position larger than index) in the value array to the next position
// move all values (with a position larger than index) in the columIndices array to the next position
Array.Copy(Values, index, Values, index + 1, ValueCount - index);
Array.Copy(ColumnIndices, index, ColumnIndices, index + 1, ValueCount - index);
Array.Copy(Values, index, Values, index + 1, valueCount - index);
Array.Copy(ColumnIndices, index, ColumnIndices, index + 1, valueCount - index);
// Add the value and the column index
Values[index] = value;
ColumnIndices[index] = column;
// increase the number of non-zero numbers by one
ValueCount += 1;
// add 1 to all the row indices for rows bigger than rowIndex
// so that they point to the correct part of the value array again.
for (var i = row + 1; i < RowPointers.Length; i++)
@ -197,10 +199,12 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
/// <remarks>WARNING: This method is not thread safe. Use "lock" with it and be sure to avoid deadlocks</remarks>
void RemoveAtIndexUnchecked(int itemIndex, int row)
{
var valueCount = RowPointers[RowPointers.Length - 1];
// Move all values (with a position larger than index) in the value array to the previous position
// move all values (with a position larger than index) in the columIndices array to the previous position
Array.Copy(Values, itemIndex + 1, Values, itemIndex, ValueCount - itemIndex - 1);
Array.Copy(ColumnIndices, itemIndex + 1, ColumnIndices, itemIndex, ValueCount - itemIndex - 1);
Array.Copy(Values, itemIndex + 1, Values, itemIndex, valueCount - itemIndex - 1);
Array.Copy(ColumnIndices, itemIndex + 1, ColumnIndices, itemIndex, valueCount - itemIndex - 1);
// Decrease value in Row
for (var i = row + 1; i < RowPointers.Length; i++)
@ -208,14 +212,14 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
RowPointers[i] -= 1;
}
ValueCount -= 1;
valueCount -= 1;
// Check whether we need to shrink the arrays. This is reasonable to do if
// there are a lot of non-zero elements and storage is two times bigger
if ((ValueCount > 1024) && (ValueCount < Values.Length / 2))
if ((valueCount > 1024) && (valueCount < Values.Length / 2))
{
Array.Resize(ref Values, ValueCount);
Array.Resize(ref ColumnIndices, ValueCount);
Array.Resize(ref Values, valueCount);
Array.Resize(ref ColumnIndices, valueCount);
}
}
@ -229,9 +233,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
public int FindItem(int row, int column)
{
// Determin bounds in columnIndices array where this item should be searched (using rowIndex)
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
return Array.BinarySearch(ColumnIndices, startIndex, endIndex - startIndex, column);
return Array.BinarySearch(ColumnIndices, RowPointers[row], RowPointers[row + 1] - RowPointers[row], column);
}
/// <summary>
@ -244,7 +246,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
int delta;
if (Values.Length > 1024)
{
delta = Values.Length / 4;
delta = Values.Length/4;
}
else
{
@ -263,7 +265,6 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
public override void Clear()
{
ValueCount = 0;
Array.Clear(RowPointers, 0, RowPointers.Length);
}
@ -275,10 +276,12 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
return;
}
var valueCount = RowPointers[RowPointers.Length - 1];
for (int row = rowIndex + rowCount - 1; row >= rowIndex; row--)
{
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
var endIndex = RowPointers[row + 1];
// empty row
if (startIndex == endIndex)
@ -297,8 +300,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
{
// Move all values (with a position larger than index) in the value array to the previous position
// move all values (with a position larger than index) in the columIndices array to the previous position
Array.Copy(Values, first + count, Values, first, ValueCount - first - count);
Array.Copy(ColumnIndices, first + count, ColumnIndices, first, ValueCount - first - count);
Array.Copy(Values, first + count, Values, first, valueCount - first - count);
Array.Copy(ColumnIndices, first + count, ColumnIndices, first, valueCount - first - count);
// Decrease value in Row
for (var k = row + 1; k < RowPointers.Length; k++)
@ -306,16 +309,16 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
RowPointers[k] -= count;
}
ValueCount -= count;
valueCount -= count;
}
}
// Check whether we need to shrink the arrays. This is reasonable to do if
// there are a lot of non-zero elements and storage is two times bigger
if ((ValueCount > 1024) && (ValueCount < Values.Length / 2))
if ((valueCount > 1024) && (valueCount < Values.Length / 2))
{
Array.Resize(ref Values, ValueCount);
Array.Resize(ref ColumnIndices, ValueCount);
Array.Resize(ref Values, valueCount);
Array.Resize(ref ColumnIndices, valueCount);
}
}
@ -350,7 +353,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
if (ValueCount != sparse.ValueCount)
{
// TODO: this is not always correct
// TODO: this is only correct if normalized
return false;
}
@ -382,7 +385,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
{
for (var i = 0; i < hashNum; i++)
{
hash = hash * 31 + values[i].GetHashCode();
hash = hash*31 + values[i].GetHashCode();
}
}
return hash;
@ -418,9 +421,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[rows] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -442,9 +445,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[rows] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -460,7 +463,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
rowPointers[row] = values.Count;
for (int col = 0; col < storage.ColumnCount; col++)
{
if (!Zero.Equals(array[row,col]))
if (!Zero.Equals(array[row, col]))
{
values.Add(array[row, col]);
columnIndices.Add(col);
@ -468,9 +471,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[storage.RowCount] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -495,9 +498,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[storage.RowCount] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -522,9 +525,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[storage.RowCount] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -551,9 +554,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[storage.RowCount] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -579,9 +582,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[storage.RowCount] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -619,9 +622,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[rows] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -654,9 +657,10 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
if (rowIterator.MoveNext()) throw new ArgumentOutOfRangeException("data", string.Format(Resources.ArgumentArrayWrongLength, rows));
}
rowPointers[rows] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -705,9 +709,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[rows] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -735,15 +739,15 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[rows] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
public static SparseCompressedRowMatrixStorage<T> OfColumnMajorList(int rows, int columns, IList<T> data)
{
if (rows * columns != data.Count)
if (rows*columns != data.Count)
{
throw new ArgumentOutOfRangeException(Resources.ArgumentMatrixDimensions);
}
@ -767,9 +771,9 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
rowPointers[rows] = values.Count;
storage.ColumnIndices = columnIndices.ToArray();
storage.Values = values.ToArray();
storage.ValueCount = values.Count;
return storage;
}
@ -803,7 +807,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
for (int row = 0; row < RowCount; row++)
{
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
var endIndex = RowPointers[row + 1];
for (var j = startIndex; j < endIndex; j++)
{
target.At(row, ColumnIndices[j], Values[j]);
@ -814,15 +818,14 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
void CopyToUnchecked(SparseCompressedRowMatrixStorage<T> target)
{
target.ValueCount = ValueCount;
target.Values = new T[ValueCount];
target.ColumnIndices = new int[ValueCount];
if (ValueCount != 0)
{
Array.Copy(Values, target.Values, ValueCount);
Buffer.BlockCopy(ColumnIndices, 0, target.ColumnIndices, 0, ValueCount * Constants.SizeOfInt);
Buffer.BlockCopy(RowPointers, 0, target.RowPointers, 0, RowCount * Constants.SizeOfInt);
Buffer.BlockCopy(ColumnIndices, 0, target.ColumnIndices, 0, ValueCount*Constants.SizeOfInt);
Buffer.BlockCopy(RowPointers, 0, target.RowPointers, 0, (RowCount + 1)*Constants.SizeOfInt);
}
}
@ -838,7 +841,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
for (int row = 0; row < RowCount; row++)
{
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
var endIndex = RowPointers[row + 1];
for (var j = startIndex; j < endIndex; j++)
{
target.At(row, ColumnIndices[j], Values[j]);
@ -860,7 +863,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
var sparseTarget = target as SparseCompressedRowMatrixStorage<T>;
if (sparseTarget != null)
{
CopySubMatrixToUnchecked(sparseTarget,
CopySubMatrixToUnchecked(sparseTarget,
sourceRowIndex, targetRowIndex, rowCount,
sourceColumnIndex, targetColumnIndex, columnCount,
skipClearing);
@ -877,7 +880,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
for (int i = sourceRowIndex, row = 0; i < sourceRowIndex + rowCount; i++, row++)
{
var startIndex = RowPointers[i];
var endIndex = i < RowPointers.Length - 1 ? RowPointers[i + 1] : ValueCount;
var endIndex = RowPointers[i + 1];
for (int j = startIndex; j < endIndex; j++)
{
@ -913,7 +916,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
rowPointers[i + rowOffset] = values.Count;
var startIndex = RowPointers[i];
var endIndex = i < RowPointers.Length - 1 ? RowPointers[i + 1] : ValueCount;
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++)
@ -927,12 +930,12 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
}
}
for(int i=targetRowIndex + rowCount; i<rowPointers.Length; i++)
for (int i = targetRowIndex + rowCount; i < rowPointers.Length; i++)
{
rowPointers[i] = values.Count;
}
target.ValueCount = values.Count;
target.RowPointers[target.RowCount] = values.Count;
target.Values = values.ToArray();
target.ColumnIndices = columnIndices.ToArray();
@ -948,7 +951,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
for (int i = sourceRowIndex, row = 0; i < sourceRowIndex + rowCount; i++, row++)
{
var startIndex = RowPointers[i];
var endIndex = i < RowPointers.Length - 1 ? RowPointers[i + 1] : ValueCount;
var endIndex = RowPointers[i + 1];
for (int j = startIndex; j < endIndex; j++)
{
@ -975,7 +978,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
// Determine bounds in columnIndices array where this item should be searched (using rowIndex)
var startIndex = RowPointers[rowIndex];
var endIndex = rowIndex < RowPointers.Length - 1 ? RowPointers[rowIndex + 1] : ValueCount;
var endIndex = RowPointers[rowIndex + 1];
if (startIndex == endIndex)
{
@ -994,14 +997,14 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
public override T[] ToRowMajorArray()
{
var ret = new T[RowCount * ColumnCount];
var ret = new T[RowCount*ColumnCount];
if (ValueCount != 0)
{
for (int row = 0; row < RowCount; row++)
{
var offset = row * ColumnCount;
var offset = row*ColumnCount;
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
var endIndex = RowPointers[row + 1];
for (var j = startIndex; j < endIndex; j++)
{
ret[offset + ColumnIndices[j]] = Values[j];
@ -1013,16 +1016,16 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
public override T[] ToColumnMajorArray()
{
var ret = new T[RowCount * ColumnCount];
var ret = new T[RowCount*ColumnCount];
if (ValueCount != 0)
{
for (int row = 0; row < RowCount; row++)
{
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
var endIndex = RowPointers[row + 1];
for (var j = startIndex; j < endIndex; j++)
{
ret[(ColumnIndices[j]) * RowCount + row] = Values[j];
ret[(ColumnIndices[j])*RowCount + row] = Values[j];
}
}
}
@ -1037,7 +1040,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
for (int row = 0; row < RowCount; row++)
{
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
var endIndex = RowPointers[row + 1];
for (var j = startIndex; j < endIndex; j++)
{
ret[row, ColumnIndices[j]] = Values[j];
@ -1056,7 +1059,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
{
for (int col = 0; col < ColumnCount; col++)
{
yield return k < (row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount) && (ColumnIndices[k]) == col
yield return k < RowPointers[row + 1] && ColumnIndices[k] == col
? Values[k++]
: Zero;
}
@ -1070,7 +1073,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
{
for (int col = 0; col < ColumnCount; col++)
{
yield return k < (row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount) && (ColumnIndices[k]) == col
yield return k < RowPointers[row + 1] && ColumnIndices[k] == col
? new Tuple<int, int, T>(row, col, Values[k++])
: new Tuple<int, int, T>(row, col, Zero);
}
@ -1087,7 +1090,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
for (int row = 0; row < RowCount; row++)
{
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
var endIndex = RowPointers[row + 1];
for (var j = startIndex; j < endIndex; j++)
{
if (!Zero.Equals(Values[j]))
@ -1102,7 +1105,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
public override void MapInplace(Func<T, T> f, bool forceMapZeros = false)
{
var newRowPointers = new int[RowCount];
var newRowPointers = new int[RowCount+1];
var newColumnIndices = new List<int>();
var newValues = new List<T>();
@ -1114,9 +1117,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
newRowPointers[row] = newValues.Count;
for (int col = 0; col < ColumnCount; col++)
{
var item = k < (row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount) && (ColumnIndices[k]) == col
? f(Values[k++])
: f(Zero);
var item = k < RowPointers[row + 1] && ColumnIndices[k] == col ? f(Values[k++]) : f(Zero);
if (!Zero.Equals(item))
{
newValues.Add(item);
@ -1131,7 +1132,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
{
newRowPointers[row] = newValues.Count;
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
var endIndex = RowPointers[row + 1];
for (var j = startIndex; j < endIndex; j++)
{
var item = f(Values[j]);
@ -1146,17 +1147,17 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
ColumnIndices = newColumnIndices.ToArray();
Values = newValues.ToArray();
ValueCount = newValues.Count;
Array.Copy(newRowPointers, RowPointers, RowCount);
newRowPointers[RowCount] = newValues.Count;
Array.Copy(newRowPointers, RowPointers, newRowPointers.Length);
}
public override void MapIndexedInplace(Func<int, int, T, T> f, bool forceMapZeros = false)
{
var newRowPointers = new int[RowCount];
var newRowPointers = new int[RowCount+1];
var newColumnIndices = new List<int>();
var newValues = new List<T>();
if (forceMapZeros || !Zero.Equals(f(0,0,Zero)))
if (forceMapZeros || !Zero.Equals(f(0, 0, Zero)))
{
int k = 0;
for (int row = 0; row < RowCount; row++)
@ -1164,9 +1165,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
newRowPointers[row] = newValues.Count;
for (int col = 0; col < ColumnCount; col++)
{
var item = k < (row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount) && (ColumnIndices[k]) == col
? f(row, col, Values[k++])
: f(row, col, Zero);
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);
@ -1181,7 +1180,7 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
{
newRowPointers[row] = newValues.Count;
var startIndex = RowPointers[row];
var endIndex = row < RowPointers.Length - 1 ? RowPointers[row + 1] : ValueCount;
var endIndex = RowPointers[row + 1];
for (var j = startIndex; j < endIndex; j++)
{
var item = f(row, ColumnIndices[j], Values[j]);
@ -1196,8 +1195,8 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
ColumnIndices = newColumnIndices.ToArray();
Values = newValues.ToArray();
ValueCount = newValues.Count;
Array.Copy(newRowPointers, RowPointers, RowCount);
newRowPointers[RowCount] = newValues.Count;
Array.Copy(newRowPointers, RowPointers, newRowPointers.Length);
}
}
}

Loading…
Cancel
Save