Browse Source

More Code Style and Cleanups

pull/729/head
Christoph Ruegg 6 years ago
parent
commit
213fd09437
  1. 102
      src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs
  2. 139
      src/Numerics/LinearAlgebra/Complex/DenseVector.cs
  3. 25
      src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs
  4. 63
      src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs
  5. 22
      src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs
  6. 42
      src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs
  7. 22
      src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs
  8. 22
      src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs
  9. 160
      src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs
  10. 253
      src/Numerics/LinearAlgebra/Complex/SparseVector.cs
  11. 115
      src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs
  12. 139
      src/Numerics/LinearAlgebra/Complex32/DenseVector.cs
  13. 25
      src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs
  14. 63
      src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs
  15. 22
      src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs
  16. 42
      src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs
  17. 22
      src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs
  18. 22
      src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs
  19. 235
      src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs
  20. 253
      src/Numerics/LinearAlgebra/Complex32/SparseVector.cs
  21. 227
      src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
  22. 134
      src/Numerics/LinearAlgebra/Double/DenseVector.cs
  23. 121
      src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs
  24. 63
      src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs
  25. 22
      src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs
  26. 42
      src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs
  27. 22
      src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs
  28. 22
      src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs
  29. 280
      src/Numerics/LinearAlgebra/Double/SparseMatrix.cs
  30. 253
      src/Numerics/LinearAlgebra/Double/SparseVector.cs
  31. 214
      src/Numerics/LinearAlgebra/Single/DenseMatrix.cs
  32. 134
      src/Numerics/LinearAlgebra/Single/DenseVector.cs
  33. 121
      src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs
  34. 63
      src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs
  35. 22
      src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs
  36. 42
      src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs
  37. 22
      src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs
  38. 22
      src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs
  39. 283
      src/Numerics/LinearAlgebra/Single/SparseMatrix.cs
  40. 253
      src/Numerics/LinearAlgebra/Single/SparseVector.cs
  41. 49
      src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs

102
src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs

@ -455,21 +455,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The matrix to store the result of the addition.</param>
protected override void DoAdd(Complex scalar, Matrix<Complex> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoAdd(scalar, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] + scalar;
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] + scalar;
}
});
base.DoAdd(scalar, result);
}
}
/// <summary>
@ -510,21 +510,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The matrix to store the result of the subtraction.</param>
protected override void DoSubtract(Complex scalar, Matrix<Complex> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoSubtract(scalar, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] - scalar;
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] - scalar;
}
});
base.DoSubtract(scalar, result);
}
}
/// <summary>
@ -563,14 +563,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The matrix to store the result of the multiplication.</param>
protected override void DoMultiply(Complex scalar, Matrix<Complex> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoMultiply(scalar, result);
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult._values);
base.DoMultiply(scalar, result);
}
}
@ -581,14 +580,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The result of the multiplication.</param>
protected override void DoMultiply(Vector<Complex> rightSide, Vector<Complex> result)
{
var denseRight = rightSide as DenseVector;
var denseResult = result as DenseVector;
if (denseRight == null || denseResult == null)
{
base.DoMultiply(rightSide, result);
}
else
if (rightSide is DenseVector denseRight && result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.MatrixMultiply(
_values,
@ -599,6 +591,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
1,
denseResult.Values);
}
else
{
base.DoMultiply(rightSide, result);
}
}
/// <summary>
@ -901,14 +897,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The matrix to store the result of the division.</param>
protected override void DoDivide(Complex divisor, Matrix<Complex> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoDivide(divisor, result);
LinearAlgebraControl.Provider.ScaleArray(1.0 / divisor, _values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(1.0/divisor, _values, denseResult._values);
base.DoDivide(divisor, result);
}
}
@ -919,16 +914,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The matrix to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Matrix<Complex> other, Matrix<Complex> result)
{
var denseOther = other as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseOther == null || denseResult == null)
if (other is DenseMatrix denseOther && result is DenseMatrix denseResult)
{
base.DoPointwiseMultiply(other, result);
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseMultiply(other, result);
}
}
@ -939,16 +931,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The matrix to store the result of the pointwise division.</param>
protected override void DoPointwiseDivide(Matrix<Complex> divisor, Matrix<Complex> result)
{
var denseDivisor = divisor as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseDivisor == null || denseResult == null)
if (divisor is DenseMatrix denseDivisor && result is DenseMatrix denseResult)
{
base.DoPointwiseDivide(divisor, result);
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseDivisor._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseDivisor._values, denseResult._values);
base.DoPointwiseDivide(divisor, result);
}
}
@ -959,16 +948,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The vector to store the result of the pointwise power.</param>
protected override void DoPointwisePower(Matrix<Complex> exponent, Matrix<Complex> result)
{
var denseExponent = exponent as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseExponent == null || denseResult == null)
if (exponent is DenseMatrix denseExponent && result is DenseMatrix denseResult)
{
base.DoPointwisePower(exponent, result);
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
base.DoPointwisePower(exponent, result);
}
}

139
src/Numerics/LinearAlgebra/Complex/DenseVector.cs

@ -206,20 +206,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The vector to store the result of the addition.</param>
protected override void DoAdd(Complex scalar, Vector<Complex> result)
{
var dense = result as DenseVector;
if (dense == null)
if (result is DenseVector dense)
{
base.DoAdd(scalar, result);
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] + scalar;
}
});
}
else
{
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] + scalar;
}
});
base.DoAdd(scalar, result);
}
}
@ -230,16 +229,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The vector to store the result of the addition.</param>
protected override void DoAdd(Vector<Complex> other, Vector<Complex> result)
{
var otherDense = other as DenseVector;
var resultDense = result as DenseVector;
if (otherDense == null || resultDense == null)
if (other is DenseVector otherDense && result is DenseVector resultDense)
{
base.DoAdd(other, result);
LinearAlgebraControl.Provider.AddArrays(_values, otherDense._values, resultDense._values);
}
else
{
LinearAlgebraControl.Provider.AddArrays(_values, otherDense._values, resultDense._values);
base.DoAdd(other, result);
}
}
@ -268,20 +264,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The vector to store the result of the subtraction.</param>
protected override void DoSubtract(Complex scalar, Vector<Complex> result)
{
var dense = result as DenseVector;
if (dense == null)
if (result is DenseVector dense)
{
base.DoSubtract(scalar, result);
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] - scalar;
}
});
}
else
{
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] - scalar;
}
});
base.DoSubtract(scalar, result);
}
}
@ -292,16 +287,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The vector to store the result of the subtraction.</param>
protected override void DoSubtract(Vector<Complex> other, Vector<Complex> result)
{
var otherDense = other as DenseVector;
var resultDense = result as DenseVector;
if (otherDense == null || resultDense == null)
if (other is DenseVector otherDense && result is DenseVector resultDense)
{
base.DoSubtract(other, result);
LinearAlgebraControl.Provider.SubtractArrays(_values, otherDense._values, resultDense._values);
}
else
{
LinearAlgebraControl.Provider.SubtractArrays(_values, otherDense._values, resultDense._values);
base.DoSubtract(other, result);
}
}
@ -345,14 +337,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">Target vector</param>
protected override void DoNegate(Vector<Complex> result)
{
var denseResult = result as DenseVector;
if (denseResult == null)
if (result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.ScaleArray(-Complex.One, _values, denseResult.Values);
}
else
{
base.DoNegate(result);
return;
}
LinearAlgebraControl.Provider.ScaleArray(-Complex.One, _values, denseResult.Values);
}
/// <summary>
@ -361,14 +353,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">Target vector</param>
protected override void DoConjugate(Vector<Complex> result)
{
var resultDense = result as DenseVector;
if (resultDense == null)
if (result is DenseVector resultDense)
{
LinearAlgebraControl.Provider.ConjugateArray(_values, resultDense._values);
}
else
{
base.DoConjugate(result);
return;
}
LinearAlgebraControl.Provider.ConjugateArray(_values, resultDense._values);
}
/// <summary>
@ -379,14 +371,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <remarks></remarks>
protected override void DoMultiply(Complex scalar, Vector<Complex> result)
{
var denseResult = result as DenseVector;
if (denseResult == null)
if (result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult.Values);
}
else
{
base.DoMultiply(scalar, result);
return;
}
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult.Values);
}
/// <summary>
@ -396,10 +388,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <returns>The sum of a[i]*b[i] for all i.</returns>
protected override Complex DoDotProduct(Vector<Complex> other)
{
var denseVector = other as DenseVector;
return denseVector == null
? base.DoDotProduct(other)
: LinearAlgebraControl.Provider.DotProduct(_values, denseVector.Values);
return other is DenseVector denseVector
? LinearAlgebraControl.Provider.DotProduct(_values, denseVector.Values)
: base.DoDotProduct(other);
}
/// <summary>
@ -409,16 +400,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <returns>The sum of conj(a[i])*b[i] for all i.</returns>
protected override Complex DoConjugateDotProduct(Vector<Complex> other)
{
var denseVector = other as DenseVector;
if (denseVector == null) return base.DoConjugateDotProduct(other);
// TODO: provide native zdotc routine
var dot = Complex.Zero;
for (var i = 0; i < _values.Length; i++)
if (other is DenseVector denseVector)
{
dot += _values[i].Conjugate()*denseVector._values[i];
var dot = Complex.Zero;
for (var i = 0; i < _values.Length; i++)
{
dot += _values[i].Conjugate() * denseVector._values[i];
}
return dot;
}
return dot;
// TODO: provide native zdotc routine
return base.DoConjugateDotProduct(other);
}
/// <summary>
@ -607,16 +601,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The vector to store the result of the pointwise division.</param>
protected override void DoPointwiseMultiply(Vector<Complex> other, Vector<Complex> result)
{
var denseOther = other as DenseVector;
var denseResult = result as DenseVector;
if (denseOther == null || denseResult == null)
if (other is DenseVector denseOther && result is DenseVector denseResult)
{
base.DoPointwiseMultiply(other, result);
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseMultiply(other, result);
}
}
@ -628,16 +619,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <remarks></remarks>
protected override void DoPointwiseDivide(Vector<Complex> divisor, Vector<Complex> result)
{
var denseDivisor = divisor as DenseVector;
var denseResult = result as DenseVector;
if (denseDivisor == null || denseResult == null)
if (divisor is DenseVector denseDivisor && result is DenseVector denseResult)
{
base.DoPointwiseDivide(divisor, result);
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseDivisor._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseDivisor._values, denseResult._values);
base.DoPointwiseDivide(divisor, result);
}
}
@ -648,16 +636,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The vector to store the result of the pointwise power.</param>
protected override void DoPointwisePower(Vector<Complex> exponent, Vector<Complex> result)
{
var denseExponent = exponent as DenseVector;
var denseResult = result as DenseVector;
if (denseExponent == null || denseResult == null)
if (exponent is DenseVector denseExponent && result is DenseVector denseResult)
{
base.DoPointwisePower(exponent, result);
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
base.DoPointwisePower(exponent, result);
}
}

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

@ -288,14 +288,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return;
}
var diagResult = result as DiagonalMatrix;
if (diagResult == null)
if (result is DiagonalMatrix diagResult)
{
base.DoMultiply(scalar, result);
LinearAlgebraControl.Provider.ScaleArray(scalar, _data, diagResult._data);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _data, diagResult._data);
base.DoMultiply(scalar, result);
}
}
@ -726,19 +725,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// this[i,i].</remarks>
public override void SetDiagonal(Vector<Complex> source)
{
var denseSource = source as DenseVector;
if (denseSource == null)
if (source is DenseVector denseSource)
{
base.SetDiagonal(source);
return;
}
if (_data.Length != denseSource.Values.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, nameof(source));
}
if (_data.Length != denseSource.Values.Length)
Array.Copy(denseSource.Values, 0, _data, 0, denseSource.Values.Length);
}
else
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, nameof(source));
base.SetDiagonal(source);
}
Array.Copy(denseSource.Values, 0, _data, 0, denseSource.Values.Length);
}
/// <summary>Calculates the induced L1 norm of this matrix.</summary>

63
src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs

@ -94,24 +94,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factor);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
var dresult = result as DenseMatrix;
if (dresult == null)
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount);
}
/// <summary>
@ -131,24 +126,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factor);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
var dresult = result as DenseVector;
if (dresult == null)
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1);
}
/// <summary>
@ -171,19 +161,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(matrix, Factor);
}
var dmatrix = matrix as DenseMatrix;
if (dmatrix == null)
if (matrix is DenseMatrix dmatrix)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
var dfactor = (DenseMatrix)Factor;
var dfactor = (DenseMatrix) Factor;
// Overwrite the existing Factor matrix with the input.
Array.Copy(dmatrix.Values, 0, dfactor.Values, 0, dmatrix.Values.Length);
// Overwrite the existing Factor matrix with the input.
Array.Copy(dmatrix.Values, 0, dfactor.Values, 0, dmatrix.Values.Length);
// Perform factorization (while overwriting).
LinearAlgebraControl.Provider.CholeskyFactor(dfactor.Values, dfactor.RowCount);
// Perform factorization (while overwriting).
LinearAlgebraControl.Provider.CholeskyFactor(dfactor.Values, dfactor.RowCount);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
}
}
}

22
src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs

@ -147,19 +147,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix)Q).Values, ((DenseMatrix)FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin);
}
/// <summary>
@ -182,19 +177,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(Q, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix)Q).Values, ((DenseMatrix)FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin);
}
}
}

42
src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs

@ -113,24 +113,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factors);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
var dresult = result as DenseMatrix;
if (dresult == null)
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
else
{
throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
/// <summary>
@ -162,24 +157,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factors);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
var dresult = result as DenseVector;
if (dresult == null)
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
else
{
throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
/// <summary>

22
src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs

@ -118,19 +118,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, Method);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, Method);
}
/// <summary>
@ -153,19 +148,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(FullR, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, Method);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, Method);
}
}
}

22
src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs

@ -105,19 +105,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values);
}
/// <summary>
@ -145,19 +140,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(VT, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values);
}
}
}

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

@ -703,51 +703,50 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
protected override void DoAdd(Matrix<Complex> other, Matrix<Complex> result)
{
var sparseOther = other as SparseMatrix;
var sparseResult = result as SparseMatrix;
if (sparseOther == null || sparseResult == null)
if (other is SparseMatrix sparseOther && result is SparseMatrix sparseResult)
{
base.DoAdd(other, result);
return;
}
if (ReferenceEquals(this, other))
{
if (!ReferenceEquals(this, result))
if (ReferenceEquals(this, other))
{
CopyTo(result);
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
LinearAlgebraControl.Provider.ScaleArray(2.0, sparseResult._storage.Values, sparseResult._storage.Values);
return;
}
LinearAlgebraControl.Provider.ScaleArray(2.0, sparseResult._storage.Values, sparseResult._storage.Values);
return;
}
SparseMatrix left;
SparseMatrix left;
if (ReferenceEquals(sparseOther, sparseResult))
{
left = this;
}
else if (ReferenceEquals(this, sparseResult))
{
left = sparseOther;
}
else
{
CopyTo(sparseResult);
left = sparseOther;
}
if (ReferenceEquals(sparseOther, sparseResult))
{
left = this;
}
else if (ReferenceEquals(this, sparseResult))
{
left = sparseOther;
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
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);
result.At(i, columnIndex, resVal);
}
}
}
else
{
CopyTo(sparseResult);
left = sparseOther;
}
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
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);
result.At(i, columnIndex, resVal);
}
base.DoAdd(other, result);
}
}
@ -835,8 +834,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return;
}
var sparseResult = result as SparseMatrix;
if (sparseResult == null)
if (result is SparseMatrix sparseResult)
{
if (!ReferenceEquals(this, result))
{
CopyTo(sparseResult);
}
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
@ -861,15 +868,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
}
}
}
else
{
if (!ReferenceEquals(this, result))
{
CopyTo(sparseResult);
}
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
}
/// <summary>
@ -1087,57 +1085,55 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">The result of the multiplication.</param>
protected override void DoTransposeAndMultiply(Matrix<Complex> other, Matrix<Complex> result)
{
var otherSparse = other as SparseMatrix;
var resultSparse = result as SparseMatrix;
if (otherSparse == null || resultSparse == null)
if (other is SparseMatrix otherSparse && result is SparseMatrix resultSparse)
{
base.DoTransposeAndMultiply(other, result);
return;
}
resultSparse.Clear();
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var otherStorage = otherSparse._storage;
resultSparse.Clear();
for (var j = 0; j < RowCount; j++)
{
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = otherStorage.RowPointers[j + 1];
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
if (startIndexOther == endIndexOther)
{
continue;
}
var otherStorage = otherSparse._storage;
for (var i = 0; i < RowCount; i++)
for (var j = 0; j < RowCount; j++)
{
// Multiply row of matrix A on row of matrix B
var startIndexThis = rowPointers[i];
var endIndexThis = rowPointers[i + 1];
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = otherStorage.RowPointers[j + 1];
if (startIndexThis == endIndexThis)
if (startIndexOther == endIndexOther)
{
continue;
}
var sum = Complex.Zero;
for (var index = startIndexOther; index < endIndexOther; index++)
for (var i = 0; i < RowCount; i++)
{
var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
if (ind >= 0)
// Multiply row of matrix A on row of matrix B
var startIndexThis = rowPointers[i];
var endIndexThis = rowPointers[i + 1];
if (startIndexThis == endIndexThis)
{
sum += otherStorage.Values[index] * values[ind];
continue;
}
var sum = Complex.Zero;
for (var index = startIndexOther; index < endIndexOther; index++)
{
var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
if (ind >= 0)
{
sum += otherStorage.Values[index] * values[ind];
}
}
}
resultSparse._storage.At(i, j, sum + result.At(i, j));
resultSparse._storage.At(i, j, sum + result.At(i, j));
}
}
}
else
{
base.DoTransposeAndMultiply(other, result);
}
}
/// <summary>

253
src/Numerics/LinearAlgebra/Complex/SparseVector.cs

@ -191,76 +191,71 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </param>
protected override void DoAdd(Vector<Complex> other, Vector<Complex> result)
{
var otherSparse = other as SparseVector;
if (otherSparse == null)
if (other is SparseVector otherSparse && result is SparseVector resultSparse)
{
base.DoAdd(other, result);
return;
}
var resultSparse = result as SparseVector;
if (resultSparse == null)
{
base.DoAdd(other, result);
return;
}
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
{
var otherValue = otherStorage.Values[j];
if (!Complex.Zero.Equals(otherValue))
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
{
var otherValue = otherStorage.Values[j];
if (!Complex.Zero.Equals(otherValue))
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], otherValue);
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], otherValue);
// TODO: result can be zero, remove?
_storage.Values[i++] += otherStorage.Values[j++];
}
else
{
i++;
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] += otherStorage.Values[j++];
}
else
{
i++;
}
}
}
else
{
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
else
{
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
{
var next = _storage.Indices[i];
if (next != last)
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
{
last = next;
result.At(next, _storage.Values[i] + otherSparse.At(next));
var next = _storage.Indices[i];
if (next != last)
{
last = next;
result.At(next, _storage.Values[i] + otherSparse.At(next));
}
i++;
}
i++;
}
else
{
var next = otherStorage.Indices[j];
if (next != last)
else
{
last = next;
result.At(next, At(next) + otherStorage.Values[j]);
var next = otherStorage.Indices[j];
if (next != last)
{
last = next;
result.At(next, At(next) + otherStorage.Values[j]);
}
j++;
}
j++;
}
}
}
else
{
base.DoAdd(other, result);
}
}
/// <summary>
@ -294,76 +289,71 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
return;
}
var otherSparse = other as SparseVector;
if (otherSparse == null)
{
base.DoSubtract(other, result);
return;
}
var resultSparse = result as SparseVector;
if (resultSparse == null)
{
base.DoSubtract(other, result);
return;
}
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
if (other is SparseVector otherSparse && result is SparseVector resultSparse)
{
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
{
var otherValue = otherStorage.Values[j];
if (!Complex.Zero.Equals(otherValue))
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
{
var otherValue = otherStorage.Values[j];
if (!Complex.Zero.Equals(otherValue))
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], -otherValue);
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], -otherValue);
// TODO: result can be zero, remove?
_storage.Values[i++] -= otherStorage.Values[j++];
}
else
{
i++;
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] -= otherStorage.Values[j++];
}
else
{
i++;
}
}
}
else
{
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
else
{
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
{
var next = _storage.Indices[i];
if (next != last)
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
{
last = next;
result.At(next, _storage.Values[i] - otherSparse.At(next));
var next = _storage.Indices[i];
if (next != last)
{
last = next;
result.At(next, _storage.Values[i] - otherSparse.At(next));
}
i++;
}
i++;
}
else
{
var next = otherStorage.Indices[j];
if (next != last)
else
{
last = next;
result.At(next, At(next) - otherStorage.Values[j]);
var next = otherStorage.Indices[j];
if (next != last)
{
last = next;
result.At(next, At(next) - otherStorage.Values[j]);
}
j++;
}
j++;
}
}
}
else
{
base.DoSubtract(other, result);
}
}
/// <summary>
@ -372,27 +362,27 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// <param name="result">Target vector</param>
protected override void DoNegate(Vector<Complex> result)
{
var sparseResult = result as SparseVector;
if (sparseResult == null)
if (result is SparseVector sparseResult)
{
if (!ReferenceEquals(this, result))
{
sparseResult._storage.ValueCount = _storage.ValueCount;
sparseResult._storage.Indices = new int[_storage.ValueCount];
Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt);
sparseResult._storage.Values = new Complex[_storage.ValueCount];
Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
}
LinearAlgebraControl.Provider.ScaleArray(-Complex.One, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], -_storage.Values[index]);
}
return;
}
if (!ReferenceEquals(this, result))
{
sparseResult._storage.ValueCount = _storage.ValueCount;
sparseResult._storage.Indices = new int[_storage.ValueCount];
Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount*Constants.SizeOfInt);
sparseResult._storage.Values = new Complex[_storage.ValueCount];
Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
}
LinearAlgebraControl.Provider.ScaleArray(-Complex.One, sparseResult._storage.Values, sparseResult._storage.Values);
}
/// <summary>
@ -434,16 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
/// </param>
protected override void DoMultiply(Complex scalar, Vector<Complex> result)
{
var sparseResult = result as SparseVector;
if (sparseResult == null)
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], scalar * _storage.Values[index]);
}
}
else
if (result is SparseVector sparseResult)
{
if (!ReferenceEquals(this, result))
{
@ -456,6 +437,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], scalar * _storage.Values[index]);
}
}
}
/// <summary>

115
src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs

@ -455,21 +455,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The matrix to store the result of the addition.</param>
protected override void DoAdd(Complex32 scalar, Matrix<Complex32> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoAdd(scalar, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] + scalar;
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] + scalar;
}
});
base.DoAdd(scalar, result);
}
}
/// <summary>
@ -510,21 +510,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The matrix to store the result of the subtraction.</param>
protected override void DoSubtract(Complex32 scalar, Matrix<Complex32> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoSubtract(scalar, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] - scalar;
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] - scalar;
}
});
base.DoSubtract(scalar, result);
}
}
/// <summary>
@ -563,14 +563,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The matrix to store the result of the multiplication.</param>
protected override void DoMultiply(Complex32 scalar, Matrix<Complex32> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoMultiply(scalar, result);
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult._values);
base.DoMultiply(scalar, result);
}
}
@ -581,14 +580,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The result of the multiplication.</param>
protected override void DoMultiply(Vector<Complex32> rightSide, Vector<Complex32> result)
{
var denseRight = rightSide as DenseVector;
var denseResult = result as DenseVector;
if (denseRight == null || denseResult == null)
{
base.DoMultiply(rightSide, result);
}
else
if (rightSide is DenseVector denseRight && result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.MatrixMultiply(
_values,
@ -599,6 +591,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
1,
denseResult.Values);
}
else
{
base.DoMultiply(rightSide, result);
}
}
/// <summary>
@ -751,14 +747,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The result of the multiplication.</param>
protected override void DoTransposeThisAndMultiply(Vector<Complex32> rightSide, Vector<Complex32> result)
{
var denseRight = rightSide as DenseVector;
var denseResult = result as DenseVector;
if (denseRight == null || denseResult == null)
{
base.DoTransposeThisAndMultiply(rightSide, result);
}
else
if (rightSide is DenseVector denseRight && result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.MatrixMultiplyWithUpdate(
Providers.LinearAlgebra.Transpose.Transpose,
@ -773,6 +762,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
0.0f,
denseResult.Values);
}
else
{
base.DoTransposeThisAndMultiply(rightSide, result);
}
}
/// <summary>
@ -905,14 +898,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The matrix to store the result of the division.</param>
protected override void DoDivide(Complex32 divisor, Matrix<Complex32> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoDivide(divisor, result);
LinearAlgebraControl.Provider.ScaleArray(1.0f / divisor, _values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(1.0f/divisor, _values, denseResult._values);
base.DoDivide(divisor, result);
}
}
@ -923,16 +915,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The matrix to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Matrix<Complex32> other, Matrix<Complex32> result)
{
var denseOther = other as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseOther == null || denseResult == null)
if (other is DenseMatrix denseOther && result is DenseMatrix denseResult)
{
base.DoPointwiseMultiply(other, result);
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseMultiply(other, result);
}
}
@ -943,16 +932,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The matrix to store the result of the pointwise division.</param>
protected override void DoPointwiseDivide(Matrix<Complex32> divisor, Matrix<Complex32> result)
{
var denseDivisor = divisor as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseDivisor == null || denseResult == null)
if (divisor is DenseMatrix denseDivisor && result is DenseMatrix denseResult)
{
base.DoPointwiseDivide(divisor, result);
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseDivisor._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseDivisor._values, denseResult._values);
base.DoPointwiseDivide(divisor, result);
}
}
@ -963,16 +949,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The vector to store the result of the pointwise power.</param>
protected override void DoPointwisePower(Matrix<Complex32> exponent, Matrix<Complex32> result)
{
var denseExponent = exponent as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseExponent == null || denseResult == null)
if (exponent is DenseMatrix denseExponent && result is DenseMatrix denseResult)
{
base.DoPointwisePower(exponent, result);
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
base.DoPointwisePower(exponent, result);
}
}

139
src/Numerics/LinearAlgebra/Complex32/DenseVector.cs

@ -206,20 +206,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The vector to store the result of the addition.</param>
protected override void DoAdd(Complex32 scalar, Vector<Complex32> result)
{
var dense = result as DenseVector;
if (dense == null)
if (result is DenseVector dense)
{
base.DoAdd(scalar, result);
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] + scalar;
}
});
}
else
{
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] + scalar;
}
});
base.DoAdd(scalar, result);
}
}
@ -230,16 +229,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The vector to store the result of the addition.</param>
protected override void DoAdd(Vector<Complex32> other, Vector<Complex32> result)
{
var otherDense = other as DenseVector;
var resultDense = result as DenseVector;
if (otherDense == null || resultDense == null)
if (other is DenseVector otherDense && result is DenseVector resultDense)
{
base.DoAdd(other, result);
LinearAlgebraControl.Provider.AddArrays(_values, otherDense._values, resultDense._values);
}
else
{
LinearAlgebraControl.Provider.AddArrays(_values, otherDense._values, resultDense._values);
base.DoAdd(other, result);
}
}
@ -268,20 +264,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The vector to store the result of the subtraction.</param>
protected override void DoSubtract(Complex32 scalar, Vector<Complex32> result)
{
var dense = result as DenseVector;
if (dense == null)
if (result is DenseVector dense)
{
base.DoSubtract(scalar, result);
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] - scalar;
}
});
}
else
{
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] - scalar;
}
});
base.DoSubtract(scalar, result);
}
}
@ -292,16 +287,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The vector to store the result of the subtraction.</param>
protected override void DoSubtract(Vector<Complex32> other, Vector<Complex32> result)
{
var otherDense = other as DenseVector;
var resultDense = result as DenseVector;
if (otherDense == null || resultDense == null)
if (other is DenseVector otherDense && result is DenseVector resultDense)
{
base.DoSubtract(other, result);
LinearAlgebraControl.Provider.SubtractArrays(_values, otherDense._values, resultDense._values);
}
else
{
LinearAlgebraControl.Provider.SubtractArrays(_values, otherDense._values, resultDense._values);
base.DoSubtract(other, result);
}
}
@ -345,14 +337,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">Target vector</param>
protected override void DoNegate(Vector<Complex32> result)
{
var denseResult = result as DenseVector;
if (denseResult == null)
if (result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.ScaleArray(-Complex32.One, _values, denseResult.Values);
}
else
{
base.DoNegate(result);
return;
}
LinearAlgebraControl.Provider.ScaleArray(-Complex32.One, _values, denseResult.Values);
}
/// <summary>
@ -361,14 +353,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">Target vector</param>
protected override void DoConjugate(Vector<Complex32> result)
{
var resultDense = result as DenseVector;
if (resultDense == null)
if (result is DenseVector resultDense)
{
LinearAlgebraControl.Provider.ConjugateArray(_values, resultDense._values);
}
else
{
base.DoConjugate(result);
return;
}
LinearAlgebraControl.Provider.ConjugateArray(_values, resultDense._values);
}
/// <summary>
@ -379,14 +371,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <remarks></remarks>
protected override void DoMultiply(Complex32 scalar, Vector<Complex32> result)
{
var denseResult = result as DenseVector;
if (denseResult == null)
if (result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult.Values);
}
else
{
base.DoMultiply(scalar, result);
return;
}
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult.Values);
}
/// <summary>
@ -396,10 +388,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <returns>The sum of a[i]*b[i] for all i.</returns>
protected override Complex32 DoDotProduct(Vector<Complex32> other)
{
var denseVector = other as DenseVector;
return denseVector == null
? base.DoDotProduct(other)
: LinearAlgebraControl.Provider.DotProduct(_values, denseVector.Values);
return other is DenseVector denseVector
? LinearAlgebraControl.Provider.DotProduct(_values, denseVector.Values)
: base.DoDotProduct(other);
}
/// <summary>
@ -409,16 +400,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <returns>The sum of conj(a[i])*b[i] for all i.</returns>
protected override Complex32 DoConjugateDotProduct(Vector<Complex32> other)
{
var denseVector = other as DenseVector;
if (denseVector == null) return base.DoConjugateDotProduct(other);
// TODO: provide native cdotc routine
var dot = Complex32.Zero;
for (var i = 0; i < _values.Length; i++)
if (other is DenseVector denseVector)
{
dot += _values[i].Conjugate() * denseVector._values[i];
var dot = Complex32.Zero;
for (var i = 0; i < _values.Length; i++)
{
dot += _values[i].Conjugate() * denseVector._values[i];
}
return dot;
}
return dot;
// TODO: provide native cdotc routine
return base.DoConjugateDotProduct(other);
}
/// <summary>
@ -607,16 +601,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The vector to store the result of the pointwise division.</param>
protected override void DoPointwiseMultiply(Vector<Complex32> other, Vector<Complex32> result)
{
var denseOther = other as DenseVector;
var denseResult = result as DenseVector;
if (denseOther == null || denseResult == null)
if (other is DenseVector denseOther && result is DenseVector denseResult)
{
base.DoPointwiseMultiply(other, result);
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseMultiply(other, result);
}
}
@ -628,16 +619,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <remarks></remarks>
protected override void DoPointwiseDivide(Vector<Complex32> divisor, Vector<Complex32> result)
{
var denseOther = divisor as DenseVector;
var denseResult = result as DenseVector;
if (denseOther == null || denseResult == null)
if (divisor is DenseVector denseOther && result is DenseVector denseResult)
{
base.DoPointwiseDivide(divisor, result);
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseDivide(divisor, result);
}
}
@ -648,16 +636,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The vector to store the result of the pointwise power.</param>
protected override void DoPointwisePower(Vector<Complex32> exponent, Vector<Complex32> result)
{
var denseExponent = exponent as DenseVector;
var denseResult = result as DenseVector;
if (denseExponent == null || denseResult == null)
if (exponent is DenseVector denseExponent && result is DenseVector denseResult)
{
base.DoPointwisePower(exponent, result);
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
base.DoPointwisePower(exponent, result);
}
}

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

@ -287,14 +287,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return;
}
var diagResult = result as DiagonalMatrix;
if (diagResult == null)
if (result is DiagonalMatrix diagResult)
{
base.DoMultiply(scalar, result);
LinearAlgebraControl.Provider.ScaleArray(scalar, _data, diagResult._data);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _data, diagResult._data);
base.DoMultiply(scalar, result);
}
}
@ -725,19 +724,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// this[i,i].</remarks>
public override void SetDiagonal(Vector<Complex32> source)
{
var denseSource = source as DenseVector;
if (denseSource == null)
if (source is DenseVector denseSource)
{
base.SetDiagonal(source);
return;
}
if (_data.Length != denseSource.Values.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, nameof(source));
}
if (_data.Length != denseSource.Values.Length)
Array.Copy(denseSource.Values, 0, _data, 0, denseSource.Values.Length);
}
else
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, nameof(source));
base.SetDiagonal(source);
}
Array.Copy(denseSource.Values, 0, _data, 0, denseSource.Values.Length);
}
/// <summary>Calculates the induced L1 norm of this matrix.</summary>

63
src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs

@ -94,24 +94,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factor);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
var dresult = result as DenseMatrix;
if (dresult == null)
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount);
}
/// <summary>
@ -131,24 +126,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factor);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
var dresult = result as DenseVector;
if (dresult == null)
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1);
}
/// <summary>
@ -171,19 +161,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(matrix, Factor);
}
var dmatrix = matrix as DenseMatrix;
if (dmatrix == null)
if (matrix is DenseMatrix dmatrix)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
var dfactor = (DenseMatrix)Factor;
var dfactor = (DenseMatrix) Factor;
// Overwrite the existing Factor matrix with the input.
Array.Copy(dmatrix.Values, 0, dfactor.Values, 0, dmatrix.Values.Length);
// Overwrite the existing Factor matrix with the input.
Array.Copy(dmatrix.Values, 0, dfactor.Values, 0, dmatrix.Values.Length);
// Perform factorization (while overwriting).
LinearAlgebraControl.Provider.CholeskyFactor(dfactor.Values, dfactor.RowCount);
// Perform factorization (while overwriting).
LinearAlgebraControl.Provider.CholeskyFactor(dfactor.Values, dfactor.RowCount);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
}
}
}

22
src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs

@ -147,19 +147,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix)Q).Values, ((DenseMatrix)FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin);
}
/// <summary>
@ -182,19 +177,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(Q, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix)Q).Values, ((DenseMatrix)FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin);
}
}
}

42
src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs

@ -113,24 +113,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factors);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
var dresult = result as DenseMatrix;
if (dresult == null)
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
else
{
throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
/// <summary>
@ -162,24 +157,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factors);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
var dresult = result as DenseVector;
if (dresult == null)
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
else
{
throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Array.Copy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length);
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
/// <summary>

22
src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs

@ -118,19 +118,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, Method);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, Method);
}
/// <summary>
@ -153,19 +148,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(FullR, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, Method);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, Method);
}
}
}

22
src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs

@ -105,19 +105,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values);
}
/// <summary>
@ -145,19 +140,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(VT, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values);
}
}
}

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

@ -703,51 +703,50 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
protected override void DoAdd(Matrix<Complex32> other, Matrix<Complex32> result)
{
var sparseOther = other as SparseMatrix;
var sparseResult = result as SparseMatrix;
if (sparseOther == null || sparseResult == null)
{
base.DoAdd(other, result);
return;
}
if (ReferenceEquals(this, other))
if (other is SparseMatrix sparseOther && result is SparseMatrix sparseResult)
{
if (!ReferenceEquals(this, result))
if (ReferenceEquals(this, other))
{
CopyTo(result);
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
LinearAlgebraControl.Provider.ScaleArray(2.0f, sparseResult._storage.Values, sparseResult._storage.Values);
return;
}
LinearAlgebraControl.Provider.ScaleArray(2.0f, sparseResult._storage.Values, sparseResult._storage.Values);
return;
}
SparseMatrix left;
SparseMatrix left;
if (ReferenceEquals(sparseOther, sparseResult))
{
left = this;
}
else if (ReferenceEquals(this, sparseResult))
{
left = sparseOther;
}
else
{
CopyTo(sparseResult);
left = sparseOther;
}
if (ReferenceEquals(sparseOther, sparseResult))
{
left = this;
}
else if (ReferenceEquals(this, sparseResult))
{
left = sparseOther;
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
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);
result.At(i, columnIndex, resVal);
}
}
}
else
{
CopyTo(sparseResult);
left = sparseOther;
}
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
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);
result.At(i, columnIndex, resVal);
}
base.DoAdd(other, result);
}
}
@ -760,59 +759,58 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
protected override void DoSubtract(Matrix<Complex32> other, Matrix<Complex32> result)
{
var sparseOther = other as SparseMatrix;
var sparseResult = result as SparseMatrix;
if (sparseOther == null || sparseResult == null)
if (other is SparseMatrix sparseOther && result is SparseMatrix sparseResult)
{
base.DoSubtract(other, result);
return;
}
if (ReferenceEquals(this, other))
{
result.Clear();
return;
}
if (ReferenceEquals(this, other))
{
result.Clear();
return;
}
var otherStorage = sparseOther._storage;
var otherStorage = sparseOther._storage;
if (ReferenceEquals(this, sparseResult))
{
for (var i = 0; i < otherStorage.RowCount; i++)
if (ReferenceEquals(this, sparseResult))
{
var endIndex = otherStorage.RowPointers[i + 1];
for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
for (var i = 0; i < otherStorage.RowCount; i++)
{
var columnIndex = otherStorage.ColumnIndices[j];
var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
result.At(i, columnIndex, resVal);
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];
result.At(i, columnIndex, resVal);
}
}
}
}
else
{
if (!ReferenceEquals(sparseOther, sparseResult))
else
{
sparseOther.CopyTo(sparseResult);
}
if (!ReferenceEquals(sparseOther, sparseResult))
{
sparseOther.CopyTo(sparseResult);
}
sparseResult.Negate(sparseResult);
sparseResult.Negate(sparseResult);
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
for (var i = 0; i < RowCount; i++)
{
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
for (var i = 0; i < RowCount; i++)
{
var columnIndex = columnIndices[j];
var resVal = sparseResult.At(i, columnIndex) + values[j];
result.At(i, columnIndex, resVal);
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];
result.At(i, columnIndex, resVal);
}
}
}
}
else
{
base.DoSubtract(other, result);
}
}
/// <summary>
@ -834,8 +832,16 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return;
}
var sparseResult = result as SparseMatrix;
if (sparseResult == null)
if (result is SparseMatrix sparseResult)
{
if (!ReferenceEquals(this, result))
{
CopyTo(sparseResult);
}
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
@ -860,15 +866,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
}
}
}
else
{
if (!ReferenceEquals(this, result))
{
CopyTo(sparseResult);
}
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
}
/// <summary>
@ -1086,57 +1083,55 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">The result of the multiplication.</param>
protected override void DoTransposeAndMultiply(Matrix<Complex32> other, Matrix<Complex32> result)
{
var otherSparse = other as SparseMatrix;
var resultSparse = result as SparseMatrix;
if (otherSparse == null || resultSparse == null)
if (other is SparseMatrix otherSparse && result is SparseMatrix resultSparse)
{
base.DoTransposeAndMultiply(other, result);
return;
}
resultSparse.Clear();
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
resultSparse.Clear();
var otherStorage = otherSparse._storage;
for (var j = 0; j < RowCount; j++)
{
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = otherStorage.RowPointers[j + 1];
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
if (startIndexOther == endIndexOther)
{
continue;
}
var otherStorage = otherSparse._storage;
for (var i = 0; i < RowCount; i++)
for (var j = 0; j < RowCount; j++)
{
// Multiply row of matrix A on row of matrix B
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = otherStorage.RowPointers[j + 1];
var startIndexThis = rowPointers[i];
var endIndexThis = rowPointers[i + 1];
if (startIndexThis == endIndexThis)
if (startIndexOther == endIndexOther)
{
continue;
}
var sum = Complex32.Zero;
for (var index = startIndexOther; index < endIndexOther; index++)
for (var i = 0; i < RowCount; i++)
{
var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
if (ind >= 0)
// Multiply row of matrix A on row of matrix B
var startIndexThis = rowPointers[i];
var endIndexThis = rowPointers[i + 1];
if (startIndexThis == endIndexThis)
{
sum += otherStorage.Values[index] * values[ind];
continue;
}
var sum = Complex32.Zero;
for (var index = startIndexOther; index < endIndexOther; index++)
{
var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
if (ind >= 0)
{
sum += otherStorage.Values[index] * values[ind];
}
}
}
resultSparse._storage.At(i, j, sum + result.At(i, j));
resultSparse._storage.At(i, j, sum + result.At(i, j));
}
}
}
else
{
base.DoTransposeAndMultiply(other, result);
}
}
/// <summary>

253
src/Numerics/LinearAlgebra/Complex32/SparseVector.cs

@ -191,76 +191,71 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </param>
protected override void DoAdd(Vector<Complex32> other, Vector<Complex32> result)
{
var otherSparse = other as SparseVector;
if (otherSparse == null)
if (other is SparseVector otherSparse && result is SparseVector resultSparse)
{
base.DoAdd(other, result);
return;
}
var resultSparse = result as SparseVector;
if (resultSparse == null)
{
base.DoAdd(other, result);
return;
}
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
{
var otherValue = otherStorage.Values[j];
if (!Complex32.Zero.Equals(otherValue))
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
{
var otherValue = otherStorage.Values[j];
if (!Complex32.Zero.Equals(otherValue))
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], otherValue);
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], otherValue);
// TODO: result can be zero, remove?
_storage.Values[i++] += otherStorage.Values[j++];
}
else
{
i++;
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] += otherStorage.Values[j++];
}
else
{
i++;
}
}
}
else
{
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
else
{
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
{
var next = _storage.Indices[i];
if (next != last)
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
{
last = next;
result.At(next, _storage.Values[i] + otherSparse.At(next));
var next = _storage.Indices[i];
if (next != last)
{
last = next;
result.At(next, _storage.Values[i] + otherSparse.At(next));
}
i++;
}
i++;
}
else
{
var next = otherStorage.Indices[j];
if (next != last)
else
{
last = next;
result.At(next, At(next) + otherStorage.Values[j]);
var next = otherStorage.Indices[j];
if (next != last)
{
last = next;
result.At(next, At(next) + otherStorage.Values[j]);
}
j++;
}
j++;
}
}
}
else
{
base.DoAdd(other, result);
}
}
/// <summary>
@ -294,76 +289,71 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
return;
}
var otherSparse = other as SparseVector;
if (otherSparse == null)
{
base.DoSubtract(other, result);
return;
}
var resultSparse = result as SparseVector;
if (resultSparse == null)
{
base.DoSubtract(other, result);
return;
}
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
if (other is SparseVector otherSparse && result is SparseVector resultSparse)
{
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
{
var otherValue = otherStorage.Values[j];
if (!Complex32.Zero.Equals(otherValue))
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
{
var otherValue = otherStorage.Values[j];
if (!Complex32.Zero.Equals(otherValue))
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], -otherValue);
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], -otherValue);
// TODO: result can be zero, remove?
_storage.Values[i++] -= otherStorage.Values[j++];
}
else
{
i++;
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] -= otherStorage.Values[j++];
}
else
{
i++;
}
}
}
else
{
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
else
{
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
{
var next = _storage.Indices[i];
if (next != last)
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
{
last = next;
result.At(next, _storage.Values[i] - otherSparse.At(next));
var next = _storage.Indices[i];
if (next != last)
{
last = next;
result.At(next, _storage.Values[i] - otherSparse.At(next));
}
i++;
}
i++;
}
else
{
var next = otherStorage.Indices[j];
if (next != last)
else
{
last = next;
result.At(next, At(next) - otherStorage.Values[j]);
var next = otherStorage.Indices[j];
if (next != last)
{
last = next;
result.At(next, At(next) - otherStorage.Values[j]);
}
j++;
}
j++;
}
}
}
else
{
base.DoSubtract(other, result);
}
}
/// <summary>
@ -372,27 +362,27 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// <param name="result">Target vector</param>
protected override void DoNegate(Vector<Complex32> result)
{
var sparseResult = result as SparseVector;
if (sparseResult == null)
if (result is SparseVector sparseResult)
{
if (!ReferenceEquals(this, result))
{
sparseResult._storage.ValueCount = _storage.ValueCount;
sparseResult._storage.Indices = new int[_storage.ValueCount];
Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt);
sparseResult._storage.Values = new Complex32[_storage.ValueCount];
Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
}
LinearAlgebraControl.Provider.ScaleArray(-Complex32.One, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], -_storage.Values[index]);
}
return;
}
if (!ReferenceEquals(this, result))
{
sparseResult._storage.ValueCount = _storage.ValueCount;
sparseResult._storage.Indices = new int[_storage.ValueCount];
Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount*Constants.SizeOfInt);
sparseResult._storage.Values = new Complex32[_storage.ValueCount];
Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
}
LinearAlgebraControl.Provider.ScaleArray(-Complex32.One, sparseResult._storage.Values, sparseResult._storage.Values);
}
/// <summary>
@ -434,16 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
/// </param>
protected override void DoMultiply(Complex32 scalar, Vector<Complex32> result)
{
var sparseResult = result as SparseVector;
if (sparseResult == null)
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], scalar * _storage.Values[index]);
}
}
else
if (result is SparseVector sparseResult)
{
if (!ReferenceEquals(this, result))
{
@ -456,6 +437,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], scalar * _storage.Values[index]);
}
}
}
/// <summary>

227
src/Numerics/LinearAlgebra/Double/DenseMatrix.cs

@ -438,21 +438,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The matrix to store the result of the addition.</param>
protected override void DoAdd(double scalar, Matrix<double> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoAdd(scalar, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] + scalar;
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] + scalar;
}
});
base.DoAdd(scalar, result);
}
}
/// <summary>
@ -493,21 +493,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The matrix to store the result of the subtraction.</param>
protected override void DoSubtract(double scalar, Matrix<double> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoSubtract(scalar, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] - scalar;
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] - scalar;
}
});
base.DoSubtract(scalar, result);
}
}
/// <summary>
@ -546,14 +546,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The matrix to store the result of the multiplication.</param>
protected override void DoMultiply(double scalar, Matrix<double> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoMultiply(scalar, result);
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult._values);
base.DoMultiply(scalar, result);
}
}
@ -564,14 +563,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The result of the multiplication.</param>
protected override void DoMultiply(Vector<double> rightSide, Vector<double> result)
{
var denseRight = rightSide as DenseVector;
var denseResult = result as DenseVector;
if (denseRight == null || denseResult == null)
{
base.DoMultiply(rightSide, result);
}
else
if (rightSide is DenseVector denseRight && result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.MatrixMultiply(
_values,
@ -582,6 +574,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
1,
denseResult.Values);
}
else
{
base.DoMultiply(rightSide, result);
}
}
/// <summary>
@ -681,14 +677,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The result of the multiplication.</param>
protected override void DoTransposeThisAndMultiply(Vector<double> rightSide, Vector<double> result)
{
var denseRight = rightSide as DenseVector;
var denseResult = result as DenseVector;
if (denseRight == null || denseResult == null)
{
base.DoTransposeThisAndMultiply(rightSide, result);
}
else
if (rightSide is DenseVector denseRight && result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.MatrixMultiplyWithUpdate(
Providers.LinearAlgebra.Transpose.Transpose,
@ -703,6 +692,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
0.0,
denseResult.Values);
}
else
{
base.DoTransposeThisAndMultiply(rightSide, result);
}
}
/// <summary>
@ -760,14 +753,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The matrix to store the result of the division.</param>
protected override void DoDivide(double divisor, Matrix<double> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoDivide(divisor, result);
LinearAlgebraControl.Provider.ScaleArray(1.0 / divisor, _values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(1.0/divisor, _values, denseResult._values);
base.DoDivide(divisor, result);
}
}
@ -778,16 +770,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The matrix to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Matrix<double> other, Matrix<double> result)
{
var denseOther = other as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseOther == null || denseResult == null)
if (other is DenseMatrix denseOther && result is DenseMatrix denseResult)
{
base.DoPointwiseMultiply(other, result);
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseMultiply(other, result);
}
}
@ -798,16 +787,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The matrix to store the result of the pointwise division.</param>
protected override void DoPointwiseDivide(Matrix<double> divisor, Matrix<double> result)
{
var denseDivisor = divisor as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseDivisor == null || denseResult == null)
if (divisor is DenseMatrix denseDivisor && result is DenseMatrix denseResult)
{
base.DoPointwiseDivide(divisor, result);
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseDivisor._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseDivisor._values, denseResult._values);
base.DoPointwiseDivide(divisor, result);
}
}
@ -818,16 +804,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The vector to store the result of the pointwise power.</param>
protected override void DoPointwisePower(Matrix<double> exponent, Matrix<double> result)
{
var denseExponent = exponent as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseExponent == null || denseResult == null)
if (exponent is DenseMatrix denseExponent && result is DenseMatrix denseResult)
{
base.DoPointwisePower(exponent, result);
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
base.DoPointwisePower(exponent, result);
}
}
@ -839,26 +822,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">Matrix to store the results in.</param>
protected override void DoModulus(double divisor, Matrix<double> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoModulus(divisor, result);
return;
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
CommonParallel.For(0, _values.Length, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = Euclid.Modulus(v[i], divisor);
}
});
}
CommonParallel.For(0, _values.Length, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = Euclid.Modulus(v[i], divisor);
}
});
base.DoModulus(divisor, result);
}
}
/// <summary>
@ -869,21 +852,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">A vector to store the results in.</param>
protected override void DoModulusByThis(double dividend, Matrix<double> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoModulusByThis(dividend, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = Euclid.Modulus(dividend, _values[i]);
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = Euclid.Modulus(dividend, _values[i]);
}
});
base.DoModulusByThis(dividend, result);
}
}
/// <summary>
@ -894,26 +877,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">Matrix to store the results in.</param>
protected override void DoRemainder(double divisor, Matrix<double> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoRemainder(divisor, result);
return;
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
CommonParallel.For(0, _values.Length, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] %= divisor;
}
});
}
CommonParallel.For(0, _values.Length, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] %= divisor;
}
});
base.DoRemainder(divisor, result);
}
}
/// <summary>
@ -924,21 +907,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">A vector to store the results in.</param>
protected override void DoRemainderByThis(double dividend, Matrix<double> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoRemainderByThis(dividend, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = dividend % _values[i];
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = dividend%_values[i];
}
});
base.DoRemainderByThis(dividend, result);
}
}
/// <summary>

134
src/Numerics/LinearAlgebra/Double/DenseVector.cs

@ -206,20 +206,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The vector to store the result of the addition.</param>
protected override void DoAdd(double scalar, Vector<double> result)
{
var dense = result as DenseVector;
if (dense == null)
if (result is DenseVector dense)
{
base.DoAdd(scalar, result);
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] + scalar;
}
});
}
else
{
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] + scalar;
}
});
base.DoAdd(scalar, result);
}
}
@ -230,16 +229,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The vector to store the result of the addition.</param>
protected override void DoAdd(Vector<double> other, Vector<double> result)
{
var otherDense = other as DenseVector;
var resultDense = result as DenseVector;
if (otherDense == null || resultDense == null)
if (other is DenseVector otherDense && result is DenseVector resultDense)
{
base.DoAdd(other, result);
LinearAlgebraControl.Provider.AddArrays(_values, otherDense._values, resultDense._values);
}
else
{
LinearAlgebraControl.Provider.AddArrays(_values, otherDense._values, resultDense._values);
base.DoAdd(other, result);
}
}
@ -278,20 +274,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The vector to store the result of the subtraction.</param>
protected override void DoSubtract(double scalar, Vector<double> result)
{
var dense = result as DenseVector;
if (dense == null)
if (result is DenseVector dense)
{
base.DoSubtract(scalar, result);
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] - scalar;
}
});
}
else
{
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] - scalar;
}
});
base.DoSubtract(scalar, result);
}
}
@ -302,16 +297,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The vector to store the result of the subtraction.</param>
protected override void DoSubtract(Vector<double> other, Vector<double> result)
{
var otherDense = other as DenseVector;
var resultDense = result as DenseVector;
if (otherDense == null || resultDense == null)
if (other is DenseVector otherDense && result is DenseVector resultDense)
{
base.DoSubtract(other, result);
LinearAlgebraControl.Provider.SubtractArrays(_values, otherDense._values, resultDense._values);
}
else
{
LinearAlgebraControl.Provider.SubtractArrays(_values, otherDense._values, resultDense._values);
base.DoSubtract(other, result);
}
}
@ -355,14 +347,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">Target vector</param>
protected override void DoNegate(Vector<double> result)
{
var denseResult = result as DenseVector;
if (denseResult == null)
if (result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.ScaleArray(-1.0d, _values, denseResult.Values);
}
else
{
base.DoNegate(result);
return;
}
LinearAlgebraControl.Provider.ScaleArray(-1.0d, _values, denseResult.Values);
}
/// <summary>
@ -373,14 +365,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <remarks></remarks>
protected override void DoMultiply(double scalar, Vector<double> result)
{
var denseResult = result as DenseVector;
if (denseResult == null)
if (result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult.Values);
}
else
{
base.DoMultiply(scalar, result);
return;
}
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult.Values);
}
/// <summary>
@ -390,10 +382,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <returns>The sum of a[i]*b[i] for all i.</returns>
protected override double DoDotProduct(Vector<double> other)
{
var denseVector = other as DenseVector;
return denseVector == null
? base.DoDotProduct(other)
: LinearAlgebraControl.Provider.DotProduct(_values, denseVector.Values);
return other is DenseVector denseVector
? LinearAlgebraControl.Provider.DotProduct(_values, denseVector.Values)
: base.DoDotProduct(other);
}
/// <summary>
@ -473,12 +464,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">A vector to store the results in.</param>
protected override void DoModulus(double divisor, Vector<double> result)
{
var dense = result as DenseVector;
if (dense == null)
{
base.DoModulus(divisor, result);
}
else
if (result is DenseVector dense)
{
CommonParallel.For(0, _length, 4096, (a, b) =>
{
@ -488,6 +474,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double
}
});
}
else
{
base.DoModulus(divisor, result);
}
}
/// <summary>
@ -498,21 +488,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">A vector to store the results in.</param>
protected override void DoRemainder(double divisor, Vector<double> result)
{
var dense = result as DenseVector;
if (dense == null)
{
base.DoRemainder(divisor, result);
}
else
if (result is DenseVector dense)
{
CommonParallel.For(0, _length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i]%divisor;
dense._values[i] = _values[i] % divisor;
}
});
}
else
{
base.DoRemainder(divisor, result);
}
}
/// <summary>
@ -689,16 +678,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The vector to store the result of the pointwise division.</param>
protected override void DoPointwiseMultiply(Vector<double> other, Vector<double> result)
{
var denseOther = other as DenseVector;
var denseResult = result as DenseVector;
if (denseOther == null || denseResult == null)
if (other is DenseVector denseOther && result is DenseVector denseResult)
{
base.DoPointwiseMultiply(other, result);
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseMultiply(other, result);
}
}
@ -710,16 +696,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <remarks></remarks>
protected override void DoPointwiseDivide(Vector<double> divisor, Vector<double> result)
{
var denseOther = divisor as DenseVector;
var denseResult = result as DenseVector;
if (denseOther == null || denseResult == null)
if (divisor is DenseVector denseOther && result is DenseVector denseResult)
{
base.DoPointwiseDivide(divisor, result);
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseDivide(divisor, result);
}
}
@ -730,16 +713,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The vector to store the result of the pointwise power.</param>
protected override void DoPointwisePower(Vector<double> exponent, Vector<double> result)
{
var denseExponent = exponent as DenseVector;
var denseResult = result as DenseVector;
if (denseExponent == null || denseResult == null)
if (exponent is DenseVector denseExponent && result is DenseVector denseResult)
{
base.DoPointwisePower(exponent, result);
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
base.DoPointwisePower(exponent, result);
}
}

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

@ -267,14 +267,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return;
}
var diagResult = result as DiagonalMatrix;
if (diagResult == null)
if (result is DiagonalMatrix diagResult)
{
base.DoMultiply(scalar, result);
LinearAlgebraControl.Provider.ScaleArray(scalar, _data, diagResult._data);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _data, diagResult._data);
base.DoMultiply(scalar, result);
}
}
@ -583,19 +582,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// this[i,i].</remarks>
public override void SetDiagonal(Vector<double> source)
{
var denseSource = source as DenseVector;
if (denseSource == null)
if (source is DenseVector denseSource)
{
base.SetDiagonal(source);
return;
}
if (_data.Length != denseSource.Values.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, nameof(source));
}
if (_data.Length != denseSource.Values.Length)
Buffer.BlockCopy(denseSource.Values, 0, _data, 0, denseSource.Values.Length * Constants.SizeOfDouble);
}
else
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, nameof(source));
base.SetDiagonal(source);
}
Buffer.BlockCopy(denseSource.Values, 0, _data, 0, denseSource.Values.Length * Constants.SizeOfDouble);
}
/// <summary>Calculates the induced L1 norm of this matrix.</summary>
@ -843,21 +842,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">Matrix to store the results in.</param>
protected override void DoModulus(double divisor, Matrix<double> result)
{
var diagonalResult = result as DiagonalMatrix;
if (diagonalResult == null)
if (result is DiagonalMatrix diagonalResult)
{
base.DoModulus(divisor, result);
return;
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = Euclid.Modulus(_data[i], divisor);
}
});
}
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
else
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = Euclid.Modulus(_data[i], divisor);
}
});
base.DoModulus(divisor, result);
}
}
/// <summary>
@ -868,21 +867,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">A vector to store the results in.</param>
protected override void DoModulusByThis(double dividend, Matrix<double> result)
{
var diagonalResult = result as DiagonalMatrix;
if (diagonalResult == null)
if (result is DiagonalMatrix diagonalResult)
{
base.DoModulusByThis(dividend, result);
return;
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = Euclid.Modulus(dividend, _data[i]);
}
});
}
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
else
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = Euclid.Modulus(dividend, _data[i]);
}
});
base.DoModulusByThis(dividend, result);
}
}
/// <summary>
@ -893,21 +892,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">Matrix to store the results in.</param>
protected override void DoRemainder(double divisor, Matrix<double> result)
{
var diagonalResult = result as DiagonalMatrix;
if (diagonalResult == null)
if (result is DiagonalMatrix diagonalResult)
{
base.DoRemainder(divisor, result);
return;
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = _data[i] % divisor;
}
});
}
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
else
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = _data[i]%divisor;
}
});
base.DoRemainder(divisor, result);
}
}
/// <summary>
@ -918,21 +917,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">A vector to store the results in.</param>
protected override void DoRemainderByThis(double dividend, Matrix<double> result)
{
var diagonalResult = result as DiagonalMatrix;
if (diagonalResult == null)
if (result is DiagonalMatrix diagonalResult)
{
base.DoRemainderByThis(dividend, result);
return;
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = dividend % _data[i];
}
});
}
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
else
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = dividend%_data[i];
}
});
base.DoRemainderByThis(dividend, result);
}
}
}
}

63
src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs

@ -92,24 +92,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factor);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfDouble);
var dresult = result as DenseMatrix;
if (dresult == null)
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfDouble);
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount);
}
/// <summary>
@ -129,24 +124,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factor);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfDouble);
var dresult = result as DenseVector;
if (dresult == null)
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfDouble);
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1);
}
/// <summary>
@ -169,19 +159,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(matrix, Factor);
}
var dmatrix = matrix as DenseMatrix;
if (dmatrix == null)
if (matrix is DenseMatrix dmatrix)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
var dfactor = (DenseMatrix)Factor;
var dfactor = (DenseMatrix) Factor;
// Overwrite the existing Factor matrix with the input.
Buffer.BlockCopy(dmatrix.Values, 0, dfactor.Values, 0, dmatrix.Values.Length * Constants.SizeOfDouble);
// Overwrite the existing Factor matrix with the input.
Buffer.BlockCopy(dmatrix.Values, 0, dfactor.Values, 0, dmatrix.Values.Length * Constants.SizeOfDouble);
// Perform factorization (while overwriting).
LinearAlgebraControl.Provider.CholeskyFactor(dfactor.Values, dfactor.RowCount);
// Perform factorization (while overwriting).
LinearAlgebraControl.Provider.CholeskyFactor(dfactor.Values, dfactor.RowCount);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
}
}
}

22
src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs

@ -145,19 +145,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix)Q).Values, ((DenseMatrix)FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin);
}
/// <summary>
@ -180,19 +175,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(Q, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix)Q).Values, ((DenseMatrix)FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin);
}
}
}

42
src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs

@ -111,24 +111,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factors);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfDouble);
var dresult = result as DenseMatrix;
if (dresult == null)
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
else
{
throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfDouble);
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
/// <summary>
@ -160,24 +155,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factors);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfDouble);
var dresult = result as DenseVector;
if (dresult == null)
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
else
{
throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfDouble);
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
/// <summary>

22
src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs

@ -116,19 +116,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, Method);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, Method);
}
/// <summary>
@ -151,19 +146,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(FullR, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, Method);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, Method);
}
}
}

22
src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs

@ -103,19 +103,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values);
}
/// <summary>
@ -143,19 +138,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(VT, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values);
}
}
}

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

@ -703,51 +703,50 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
protected override void DoAdd(Matrix<double> other, Matrix<double> result)
{
var sparseOther = other as SparseMatrix;
var sparseResult = result as SparseMatrix;
if (sparseOther == null || sparseResult == null)
{
base.DoAdd(other, result);
return;
}
if (ReferenceEquals(this, other))
if (other is SparseMatrix sparseOther && result is SparseMatrix sparseResult)
{
if (!ReferenceEquals(this, result))
if (ReferenceEquals(this, other))
{
CopyTo(result);
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
LinearAlgebraControl.Provider.ScaleArray(2.0, sparseResult._storage.Values, sparseResult._storage.Values);
return;
}
LinearAlgebraControl.Provider.ScaleArray(2.0, sparseResult._storage.Values, sparseResult._storage.Values);
return;
}
SparseMatrix left;
SparseMatrix left;
if (ReferenceEquals(sparseOther, sparseResult))
{
left = this;
}
else if (ReferenceEquals(this, sparseResult))
{
left = sparseOther;
}
else
{
CopyTo(sparseResult);
left = sparseOther;
}
if (ReferenceEquals(sparseOther, sparseResult))
{
left = this;
}
else if (ReferenceEquals(this, sparseResult))
{
left = sparseOther;
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
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);
result.At(i, columnIndex, resVal);
}
}
}
else
{
CopyTo(sparseResult);
left = sparseOther;
}
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
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);
result.At(i, columnIndex, resVal);
}
base.DoAdd(other, result);
}
}
@ -760,60 +759,58 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
protected override void DoSubtract(Matrix<double> other, Matrix<double> result)
{
var sparseOther = other as SparseMatrix;
var sparseResult = result as SparseMatrix;
if (sparseOther == null || sparseResult == null)
{
base.DoSubtract(other, result);
return;
}
if (ReferenceEquals(this, other))
if (other is SparseMatrix sparseOther && result is SparseMatrix sparseResult)
{
result.Clear();
return;
}
if (ReferenceEquals(this, other))
{
result.Clear();
return;
}
var otherStorage = sparseOther._storage;
var otherStorage = sparseOther._storage;
if (ReferenceEquals(this, sparseResult))
{
for (var i = 0; i < otherStorage.RowCount; i++)
if (ReferenceEquals(this, sparseResult))
{
var endIndex = otherStorage.RowPointers[i + 1];
for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
for (var i = 0; i < otherStorage.RowCount; i++)
{
var columnIndex = otherStorage.ColumnIndices[j];
var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
result.At(i, columnIndex, resVal);
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];
result.At(i, columnIndex, resVal);
}
}
}
}
else
{
if (!ReferenceEquals(sparseOther, sparseResult))
else
{
sparseOther.CopyTo(sparseResult);
}
if (!ReferenceEquals(sparseOther, sparseResult))
{
sparseOther.CopyTo(sparseResult);
}
sparseResult.Negate(sparseResult);
sparseResult.Negate(sparseResult);
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
for (var i = 0; i < RowCount; i++)
{
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
for (var i = 0; i < RowCount; i++)
{
var columnIndex = columnIndices[j];
var resVal = sparseResult.At(i, columnIndex) + values[j];
result.At(i, columnIndex, resVal);
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];
result.At(i, columnIndex, resVal);
}
}
}
}
else
{
base.DoSubtract(other, result);
}
}
/// <summary>
@ -835,8 +832,16 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return;
}
var sparseResult = result as SparseMatrix;
if (sparseResult == null)
if (result is SparseMatrix sparseResult)
{
if (!ReferenceEquals(this, result))
{
CopyTo(sparseResult);
}
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
@ -861,15 +866,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
}
}
}
else
{
if (!ReferenceEquals(this, result))
{
CopyTo(sparseResult);
}
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
}
/// <summary>
@ -1086,55 +1082,53 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">The result of the multiplication.</param>
protected override void DoTransposeAndMultiply(Matrix<double> other, Matrix<double> result)
{
var otherSparse = other as SparseMatrix;
var resultSparse = result as SparseMatrix;
if (otherSparse == null || resultSparse == null)
if (other is SparseMatrix otherSparse && result is SparseMatrix resultSparse)
{
base.DoTransposeAndMultiply(other, result);
return;
}
resultSparse.Clear();
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
resultSparse.Clear();
var otherStorage = otherSparse._storage;
for (var j = 0; j < RowCount; j++)
{
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = otherStorage.RowPointers[j + 1];
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
if (startIndexOther == endIndexOther)
{
continue;
}
var otherStorage = otherSparse._storage;
for (var i = 0; i < RowCount; i++)
for (var j = 0; j < RowCount; j++)
{
var startIndexThis = rowPointers[i];
var endIndexThis = rowPointers[i + 1];
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = otherStorage.RowPointers[j + 1];
if (startIndexThis == endIndexThis)
if (startIndexOther == endIndexOther)
{
continue;
}
var sum = 0d;
for (var index = startIndexOther; index < endIndexOther; index++)
for (var i = 0; i < RowCount; i++)
{
var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
if (ind >= 0)
var startIndexThis = rowPointers[i];
var endIndexThis = rowPointers[i + 1];
if (startIndexThis == endIndexThis)
{
sum += otherStorage.Values[index]*values[ind];
continue;
}
}
resultSparse._storage.At(i, j, sum + result.At(i, j));
var sum = 0d;
for (var index = startIndexOther; index < endIndexOther; index++)
{
var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
if (ind >= 0)
{
sum += otherStorage.Values[index] * values[ind];
}
}
resultSparse._storage.At(i, j, sum + result.At(i, j));
}
}
}
else
{
base.DoTransposeAndMultiply(other, result);
}
}
/// <summary>
@ -1261,22 +1255,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">Matrix to store the results in.</param>
protected override void DoModulus(double divisor, Matrix<double> result)
{
var sparseResult = result as SparseMatrix;
if (sparseResult == null)
if (result is SparseMatrix sparseResult)
{
base.DoModulus(divisor, result);
return;
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
var resultStorage = sparseResult._storage;
for (var index = 0; index < resultStorage.Values.Length; index++)
{
resultStorage.Values[index] = Euclid.Modulus(resultStorage.Values[index], divisor);
}
}
var resultStorage = sparseResult._storage;
for (var index = 0; index < resultStorage.Values.Length; index++)
else
{
resultStorage.Values[index] = Euclid.Modulus(resultStorage.Values[index], divisor);
base.DoModulus(divisor, result);
}
}
@ -1288,22 +1282,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">Matrix to store the results in.</param>
protected override void DoRemainder(double divisor, Matrix<double> result)
{
var sparseResult = result as SparseMatrix;
if (sparseResult == null)
if (result is SparseMatrix sparseResult)
{
base.DoRemainder(divisor, result);
return;
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
var resultStorage = sparseResult._storage;
for (var index = 0; index < resultStorage.Values.Length; index++)
{
resultStorage.Values[index] %= divisor;
}
}
var resultStorage = sparseResult._storage;
for (var index = 0; index < resultStorage.Values.Length; index++)
else
{
resultStorage.Values[index] %= divisor;
base.DoRemainder(divisor, result);
}
}

253
src/Numerics/LinearAlgebra/Double/SparseVector.cs

@ -191,76 +191,71 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </param>
protected override void DoAdd(Vector<double> other, Vector<double> result)
{
var otherSparse = other as SparseVector;
if (otherSparse == null)
if (other is SparseVector otherSparse && result is SparseVector resultSparse)
{
base.DoAdd(other, result);
return;
}
var resultSparse = result as SparseVector;
if (resultSparse == null)
{
base.DoAdd(other, result);
return;
}
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
{
var otherValue = otherStorage.Values[j];
if (otherValue != 0.0)
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
{
var otherValue = otherStorage.Values[j];
if (otherValue != 0.0)
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], otherValue);
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], otherValue);
// TODO: result can be zero, remove?
_storage.Values[i++] += otherStorage.Values[j++];
}
else
{
i++;
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] += otherStorage.Values[j++];
}
else
{
i++;
}
}
}
else
{
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
else
{
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
{
var next = _storage.Indices[i];
if (next != last)
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
{
last = next;
result.At(next, _storage.Values[i] + otherSparse.At(next));
var next = _storage.Indices[i];
if (next != last)
{
last = next;
result.At(next, _storage.Values[i] + otherSparse.At(next));
}
i++;
}
i++;
}
else
{
var next = otherStorage.Indices[j];
if (next != last)
else
{
last = next;
result.At(next, At(next) + otherStorage.Values[j]);
var next = otherStorage.Indices[j];
if (next != last)
{
last = next;
result.At(next, At(next) + otherStorage.Values[j]);
}
j++;
}
j++;
}
}
}
else
{
base.DoAdd(other, result);
}
}
/// <summary>
@ -294,76 +289,71 @@ namespace MathNet.Numerics.LinearAlgebra.Double
return;
}
var otherSparse = other as SparseVector;
if (otherSparse == null)
if (other is SparseVector otherSparse && result is SparseVector resultSparse)
{
base.DoSubtract(other, result);
return;
}
var resultSparse = result as SparseVector;
if (resultSparse == null)
{
base.DoSubtract(other, result);
return;
}
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
{
var otherValue = otherStorage.Values[j];
if (otherValue != 0.0)
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], -otherValue);
var otherValue = otherStorage.Values[j];
if (otherValue != 0.0)
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], -otherValue);
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] -= otherStorage.Values[j++];
}
else
{
i++;
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] -= otherStorage.Values[j++];
}
else
{
i++;
}
}
}
else
{
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
else
{
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
{
var next = _storage.Indices[i];
if (next != last)
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
{
last = next;
result.At(next, _storage.Values[i] - otherSparse.At(next));
var next = _storage.Indices[i];
if (next != last)
{
last = next;
result.At(next, _storage.Values[i] - otherSparse.At(next));
}
i++;
}
i++;
}
else
{
var next = otherStorage.Indices[j];
if (next != last)
else
{
last = next;
result.At(next, At(next) - otherStorage.Values[j]);
var next = otherStorage.Indices[j];
if (next != last)
{
last = next;
result.At(next, At(next) - otherStorage.Values[j]);
}
j++;
}
j++;
}
}
}
else
{
base.DoSubtract(other, result);
}
}
/// <summary>
@ -372,27 +362,27 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// <param name="result">Target vector</param>
protected override void DoNegate(Vector<double> result)
{
var sparseResult = result as SparseVector;
if (sparseResult == null)
if (result is SparseVector sparseResult)
{
if (!ReferenceEquals(this, result))
{
sparseResult._storage.ValueCount = _storage.ValueCount;
sparseResult._storage.Indices = new int[_storage.ValueCount];
Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt);
sparseResult._storage.Values = new double[_storage.ValueCount];
Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
}
LinearAlgebraControl.Provider.ScaleArray(-1.0d, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], -_storage.Values[index]);
}
return;
}
if (!ReferenceEquals(this, result))
{
sparseResult._storage.ValueCount = _storage.ValueCount;
sparseResult._storage.Indices = new int[_storage.ValueCount];
Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt);
sparseResult._storage.Values = new double[_storage.ValueCount];
Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
}
LinearAlgebraControl.Provider.ScaleArray(-1.0d, sparseResult._storage.Values, sparseResult._storage.Values);
}
/// <summary>
@ -406,16 +396,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double
/// </param>
protected override void DoMultiply(double scalar, Vector<double> result)
{
var sparseResult = result as SparseVector;
if (sparseResult == null)
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], scalar * _storage.Values[index]);
}
}
else
if (result is SparseVector sparseResult)
{
if (!ReferenceEquals(this, result))
{
@ -428,6 +409,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], scalar * _storage.Values[index]);
}
}
}
/// <summary>

214
src/Numerics/LinearAlgebra/Single/DenseMatrix.cs

@ -438,21 +438,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The matrix to store the result of the addition.</param>
protected override void DoAdd(float scalar, Matrix<float> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoAdd(scalar, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] + scalar;
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] + scalar;
}
});
base.DoAdd(scalar, result);
}
}
/// <summary>
@ -493,21 +493,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The matrix to store the result of the subtraction.</param>
protected override void DoSubtract(float scalar, Matrix<float> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoSubtract(scalar, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] - scalar;
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = _values[i] - scalar;
}
});
base.DoSubtract(scalar, result);
}
}
/// <summary>
@ -546,14 +546,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The matrix to store the result of the multiplication.</param>
protected override void DoMultiply(float scalar, Matrix<float> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoMultiply(scalar, result);
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult._values);
base.DoMultiply(scalar, result);
}
}
@ -564,14 +563,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The result of the multiplication.</param>
protected override void DoMultiply(Vector<float> rightSide, Vector<float> result)
{
var denseRight = rightSide as DenseVector;
var denseResult = result as DenseVector;
if (denseRight == null || denseResult == null)
{
base.DoMultiply(rightSide, result);
}
else
if (rightSide is DenseVector denseRight && result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.MatrixMultiply(
_values,
@ -582,6 +574,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
1,
denseResult.Values);
}
else
{
base.DoMultiply(rightSide, result);
}
}
/// <summary>
@ -757,14 +753,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The matrix to store the result of the division.</param>
protected override void DoDivide(float divisor, Matrix<float> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoDivide(divisor, result);
LinearAlgebraControl.Provider.ScaleArray(1.0f / divisor, _values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(1.0f/divisor, _values, denseResult._values);
base.DoDivide(divisor, result);
}
}
@ -775,16 +770,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The matrix to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Matrix<float> other, Matrix<float> result)
{
var denseOther = other as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseOther == null || denseResult == null)
if (other is DenseMatrix denseOther && result is DenseMatrix denseResult)
{
base.DoPointwiseMultiply(other, result);
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseMultiply(other, result);
}
}
@ -795,16 +787,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The matrix to store the result of the pointwise division.</param>
protected override void DoPointwiseDivide(Matrix<float> divisor, Matrix<float> result)
{
var denseOther = divisor as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseOther == null || denseResult == null)
if (divisor is DenseMatrix denseOther && result is DenseMatrix denseResult)
{
base.DoPointwiseDivide(divisor, result);
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseDivide(divisor, result);
}
}
@ -815,16 +804,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The vector to store the result of the pointwise power.</param>
protected override void DoPointwisePower(Matrix<float> exponent, Matrix<float> result)
{
var denseExponent = exponent as DenseMatrix;
var denseResult = result as DenseMatrix;
if (denseExponent == null || denseResult == null)
if (exponent is DenseMatrix denseExponent && result is DenseMatrix denseResult)
{
base.DoPointwisePower(exponent, result);
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
base.DoPointwisePower(exponent, result);
}
}
@ -836,26 +822,26 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">Matrix to store the results in.</param>
protected override void DoModulus(float divisor, Matrix<float> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoModulus(divisor, result);
return;
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
CommonParallel.For(0, _values.Length, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = Euclid.Modulus(v[i], divisor);
}
});
}
CommonParallel.For(0, _values.Length, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = Euclid.Modulus(v[i], divisor);
}
});
base.DoModulus(divisor, result);
}
}
/// <summary>
@ -866,21 +852,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">A vector to store the results in.</param>
protected override void DoModulusByThis(float dividend, Matrix<float> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoModulusByThis(dividend, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = Euclid.Modulus(dividend, _values[i]);
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = Euclid.Modulus(dividend, _values[i]);
}
});
base.DoModulusByThis(dividend, result);
}
}
/// <summary>
@ -891,26 +877,26 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">Matrix to store the results in.</param>
protected override void DoRemainder(float divisor, Matrix<float> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoRemainder(divisor, result);
return;
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
CommonParallel.For(0, _values.Length, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] %= divisor;
}
});
}
CommonParallel.For(0, _values.Length, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] %= divisor;
}
});
base.DoRemainder(divisor, result);
}
}
/// <summary>
@ -921,21 +907,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">A vector to store the results in.</param>
protected override void DoRemainderByThis(float dividend, Matrix<float> result)
{
var denseResult = result as DenseMatrix;
if (denseResult == null)
if (result is DenseMatrix denseResult)
{
base.DoRemainderByThis(dividend, result);
return;
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = dividend % _values[i];
}
});
}
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
else
{
var v = denseResult._values;
for (int i = a; i < b; i++)
{
v[i] = dividend%_values[i];
}
});
base.DoRemainderByThis(dividend, result);
}
}
/// <summary>

134
src/Numerics/LinearAlgebra/Single/DenseVector.cs

@ -206,20 +206,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The vector to store the result of the addition.</param>
protected override void DoAdd(float scalar, Vector<float> result)
{
var dense = result as DenseVector;
if (dense == null)
if (result is DenseVector dense)
{
base.DoAdd(scalar, result);
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] + scalar;
}
});
}
else
{
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] + scalar;
}
});
base.DoAdd(scalar, result);
}
}
@ -230,16 +229,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The vector to store the result of the addition.</param>
protected override void DoAdd(Vector<float> other, Vector<float> result)
{
var otherDense = other as DenseVector;
var resultDense = result as DenseVector;
if (otherDense == null || resultDense == null)
if (other is DenseVector otherDense && result is DenseVector resultDense)
{
base.DoAdd(other, result);
LinearAlgebraControl.Provider.AddArrays(_values, otherDense._values, resultDense._values);
}
else
{
LinearAlgebraControl.Provider.AddArrays(_values, otherDense._values, resultDense._values);
base.DoAdd(other, result);
}
}
@ -268,20 +264,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The vector to store the result of the subtraction.</param>
protected override void DoSubtract(float scalar, Vector<float> result)
{
var dense = result as DenseVector;
if (dense == null)
if (result is DenseVector dense)
{
base.DoSubtract(scalar, result);
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] - scalar;
}
});
}
else
{
CommonParallel.For(0, _values.Length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i] - scalar;
}
});
base.DoSubtract(scalar, result);
}
}
@ -292,16 +287,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The vector to store the result of the subtraction.</param>
protected override void DoSubtract(Vector<float> other, Vector<float> result)
{
var otherDense = other as DenseVector;
var resultDense = result as DenseVector;
if (otherDense == null || resultDense == null)
if (other is DenseVector otherDense && result is DenseVector resultDense)
{
base.DoSubtract(other, result);
LinearAlgebraControl.Provider.SubtractArrays(_values, otherDense._values, resultDense._values);
}
else
{
LinearAlgebraControl.Provider.SubtractArrays(_values, otherDense._values, resultDense._values);
base.DoSubtract(other, result);
}
}
@ -345,14 +337,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">Target vector</param>
protected override void DoNegate(Vector<float> result)
{
var denseResult = result as DenseVector;
if (denseResult == null)
if (result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.ScaleArray(-1.0f, _values, denseResult.Values);
}
else
{
base.DoNegate(result);
return;
}
LinearAlgebraControl.Provider.ScaleArray(-1.0f, _values, denseResult.Values);
}
/// <summary>
@ -363,14 +355,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <remarks></remarks>
protected override void DoMultiply(float scalar, Vector<float> result)
{
var denseResult = result as DenseVector;
if (denseResult == null)
if (result is DenseVector denseResult)
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult.Values);
}
else
{
base.DoMultiply(scalar, result);
return;
}
LinearAlgebraControl.Provider.ScaleArray(scalar, _values, denseResult.Values);
}
/// <summary>
@ -380,10 +372,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <returns>The sum of a[i]*b[i] for all i.</returns>
protected override float DoDotProduct(Vector<float> other)
{
var denseVector = other as DenseVector;
return denseVector == null
? base.DoDotProduct(other)
: LinearAlgebraControl.Provider.DotProduct(_values, denseVector.Values);
return other is DenseVector denseVector
? LinearAlgebraControl.Provider.DotProduct(_values, denseVector.Values)
: base.DoDotProduct(other);
}
/// <summary>
@ -463,12 +454,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">A vector to store the results in.</param>
protected override void DoModulus(float divisor, Vector<float> result)
{
var dense = result as DenseVector;
if (dense == null)
{
base.DoModulus(divisor, result);
}
else
if (result is DenseVector dense)
{
CommonParallel.For(0, _length, 4096, (a, b) =>
{
@ -478,6 +464,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single
}
});
}
else
{
base.DoModulus(divisor, result);
}
}
/// <summary>
@ -488,21 +478,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">A vector to store the results in.</param>
protected override void DoRemainder(float divisor, Vector<float> result)
{
var dense = result as DenseVector;
if (dense == null)
{
base.DoRemainder(divisor, result);
}
else
if (result is DenseVector dense)
{
CommonParallel.For(0, _length, 4096, (a, b) =>
{
for (int i = a; i < b; i++)
{
dense._values[i] = _values[i]%divisor;
dense._values[i] = _values[i] % divisor;
}
});
}
else
{
base.DoRemainder(divisor, result);
}
}
/// <summary>
@ -680,16 +669,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The vector to store the result of the pointwise multiplication.</param>
protected override void DoPointwiseMultiply(Vector<float> other, Vector<float> result)
{
var denseOther = other as DenseVector;
var denseResult = result as DenseVector;
if (denseOther == null || denseResult == null)
if (other is DenseVector denseOther && result is DenseVector denseResult)
{
base.DoPointwiseMultiply(other, result);
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseMultiplyArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseMultiply(other, result);
}
}
@ -701,16 +687,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <remarks></remarks>
protected override void DoPointwiseDivide(Vector<float> divisor, Vector<float> result)
{
var denseOther = divisor as DenseVector;
var denseResult = result as DenseVector;
if (denseOther == null || denseResult == null)
if (divisor is DenseVector denseOther && result is DenseVector denseResult)
{
base.DoPointwiseDivide(divisor, result);
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseOther._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWiseDivideArrays(_values, denseOther._values, denseResult._values);
base.DoPointwiseDivide(divisor, result);
}
}
@ -721,16 +704,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The vector to store the result of the pointwise power.</param>
protected override void DoPointwisePower(Vector<float> exponent, Vector<float> result)
{
var denseExponent = exponent as DenseVector;
var denseResult = result as DenseVector;
if (denseExponent == null || denseResult == null)
if (exponent is DenseVector denseExponent && result is DenseVector denseResult)
{
base.DoPointwisePower(exponent, result);
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
}
else
{
LinearAlgebraControl.Provider.PointWisePowerArrays(_values, denseExponent._values, denseResult._values);
base.DoPointwisePower(exponent, result);
}
}

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

@ -267,14 +267,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return;
}
var diagResult = result as DiagonalMatrix;
if (diagResult == null)
if (result is DiagonalMatrix diagResult)
{
base.DoMultiply(scalar, result);
LinearAlgebraControl.Provider.ScaleArray(scalar, _data, diagResult._data);
}
else
{
LinearAlgebraControl.Provider.ScaleArray(scalar, _data, diagResult._data);
base.DoMultiply(scalar, result);
}
}
@ -583,19 +582,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// this[i,i].</remarks>
public override void SetDiagonal(Vector<float> source)
{
var denseSource = source as DenseVector;
if (denseSource == null)
if (source is DenseVector denseSource)
{
base.SetDiagonal(source);
return;
}
if (_data.Length != denseSource.Values.Length)
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, nameof(source));
}
if (_data.Length != denseSource.Values.Length)
Buffer.BlockCopy(denseSource.Values, 0, _data, 0, denseSource.Values.Length * Constants.SizeOfFloat);
}
else
{
throw new ArgumentException(Resources.ArgumentVectorsSameLength, nameof(source));
base.SetDiagonal(source);
}
Buffer.BlockCopy(denseSource.Values, 0, _data, 0, denseSource.Values.Length * Constants.SizeOfFloat);
}
/// <summary>Calculates the induced L1 norm of this matrix.</summary>
@ -843,21 +842,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">Matrix to store the results in.</param>
protected override void DoModulus(float divisor, Matrix<float> result)
{
var diagonalResult = result as DiagonalMatrix;
if (diagonalResult == null)
if (result is DiagonalMatrix diagonalResult)
{
base.DoModulus(divisor, result);
return;
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = Euclid.Modulus(_data[i], divisor);
}
});
}
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
else
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = Euclid.Modulus(_data[i], divisor);
}
});
base.DoModulus(divisor, result);
}
}
/// <summary>
@ -868,21 +867,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">A vector to store the results in.</param>
protected override void DoModulusByThis(float dividend, Matrix<float> result)
{
var diagonalResult = result as DiagonalMatrix;
if (diagonalResult == null)
if (result is DiagonalMatrix diagonalResult)
{
base.DoModulusByThis(dividend, result);
return;
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = Euclid.Modulus(dividend, _data[i]);
}
});
}
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
else
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = Euclid.Modulus(dividend, _data[i]);
}
});
base.DoModulusByThis(dividend, result);
}
}
/// <summary>
@ -893,21 +892,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">Matrix to store the results in.</param>
protected override void DoRemainder(float divisor, Matrix<float> result)
{
var diagonalResult = result as DiagonalMatrix;
if (diagonalResult == null)
if (result is DiagonalMatrix diagonalResult)
{
base.DoRemainder(divisor, result);
return;
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = _data[i] % divisor;
}
});
}
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
else
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = _data[i]%divisor;
}
});
base.DoRemainder(divisor, result);
}
}
/// <summary>
@ -918,21 +917,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">A vector to store the results in.</param>
protected override void DoRemainderByThis(float dividend, Matrix<float> result)
{
var diagonalResult = result as DiagonalMatrix;
if (diagonalResult == null)
if (result is DiagonalMatrix diagonalResult)
{
base.DoRemainderByThis(dividend, result);
return;
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = dividend % _data[i];
}
});
}
CommonParallel.For(0, _data.Length, 4096, (a, b) =>
else
{
var r = diagonalResult._data;
for (var i = a; i < b; i++)
{
r[i] = dividend%_data[i];
}
});
base.DoRemainderByThis(dividend, result);
}
}
}
}

63
src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs

@ -92,24 +92,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factor);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfFloat);
var dresult = result as DenseMatrix;
if (dresult == null)
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfFloat);
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, dresult.ColumnCount);
}
/// <summary>
@ -129,24 +124,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factor);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfFloat);
var dresult = result as DenseVector;
if (dresult == null)
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfFloat);
// Cholesky solve by overwriting result.
var dfactor = (DenseMatrix) Factor;
LinearAlgebraControl.Provider.CholeskySolveFactored(dfactor.Values, dfactor.RowCount, dresult.Values, 1);
}
/// <summary>
@ -169,19 +159,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(matrix, Factor);
}
var dmatrix = matrix as DenseMatrix;
if (dmatrix == null)
if (matrix is DenseMatrix dmatrix)
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
var dfactor = (DenseMatrix)Factor;
var dfactor = (DenseMatrix) Factor;
// Overwrite the existing Factor matrix with the input.
Buffer.BlockCopy(dmatrix.Values, 0, dfactor.Values, 0, dmatrix.Values.Length * Constants.SizeOfFloat);
// Overwrite the existing Factor matrix with the input.
Buffer.BlockCopy(dmatrix.Values, 0, dfactor.Values, 0, dmatrix.Values.Length * Constants.SizeOfFloat);
// Perform factorization (while overwriting).
LinearAlgebraControl.Provider.CholeskyFactor(dfactor.Values, dfactor.RowCount);
// Perform factorization (while overwriting).
LinearAlgebraControl.Provider.CholeskyFactor(dfactor.Values, dfactor.RowCount);
}
else
{
throw new NotSupportedException("Can only do Cholesky factorization for dense matrices at the moment.");
}
}
}
}

22
src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs

@ -145,19 +145,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix)Q).Values, ((DenseMatrix)FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, input.ColumnCount, dresult.Values, QRMethod.Thin);
}
/// <summary>
@ -180,19 +175,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(Q, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do GramSchmidt factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix)Q).Values, ((DenseMatrix)FullR).Values, Q.RowCount, FullR.ColumnCount, null, dinput.Values, 1, dresult.Values, QRMethod.Thin);
}
}
}

42
src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs

@ -111,24 +111,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factors);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfFloat);
var dresult = result as DenseMatrix;
if (dresult == null)
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
else
{
throw new NotSupportedException("Can only do LU factorization for dense matrices at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfFloat);
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(input.ColumnCount, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
/// <summary>
@ -160,24 +155,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(input, Factors);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length * Constants.SizeOfFloat);
var dresult = result as DenseVector;
if (dresult == null)
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
else
{
throw new NotSupportedException("Can only do LU factorization for dense vectors at the moment.");
}
// Copy the contents of input to result.
Buffer.BlockCopy(dinput.Values, 0, dresult.Values, 0, dinput.Values.Length*Constants.SizeOfFloat);
// LU solve by overwriting result.
var dfactors = (DenseMatrix) Factors;
LinearAlgebraControl.Provider.LUSolveFactored(1, dfactors.Values, dfactors.RowCount, Pivots, dresult.Values);
}
/// <summary>

22
src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs

@ -116,19 +116,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, Method);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do QR factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, input.ColumnCount, dresult.Values, Method);
}
/// <summary>
@ -151,19 +146,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(FullR, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, Method);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do QR factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.QRSolveFactored(((DenseMatrix) Q).Values, ((DenseMatrix) FullR).Values, Q.RowCount, FullR.ColumnCount, Tau, dinput.Values, 1, dresult.Values, Method);
}
}
}

22
src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs

@ -104,19 +104,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw new ArgumentException(Resources.ArgumentMatrixSameColumnDimension);
}
var dinput = input as DenseMatrix;
if (dinput == null)
if (input is DenseMatrix dinput && result is DenseMatrix dresult)
{
throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment.");
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values);
}
var dresult = result as DenseMatrix;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do SVD factorization for dense matrices at the moment.");
}
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, input.ColumnCount, dresult.Values);
}
/// <summary>
@ -144,19 +139,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
throw Matrix.DimensionsDontMatch<ArgumentException>(VT, result);
}
var dinput = input as DenseVector;
if (dinput == null)
if (input is DenseVector dinput && result is DenseVector dresult)
{
throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment.");
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values);
}
var dresult = result as DenseVector;
if (dresult == null)
else
{
throw new NotSupportedException("Can only do SVD factorization for dense vectors at the moment.");
}
LinearAlgebraControl.Provider.SvdSolveFactored(U.RowCount, VT.ColumnCount, ((DenseVector) S).Values, ((DenseMatrix) U).Values, ((DenseMatrix) VT).Values, dinput.Values, 1, dresult.Values);
}
}
}

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

@ -708,51 +708,50 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
protected override void DoAdd(Matrix<float> other, Matrix<float> result)
{
var sparseOther = other as SparseMatrix;
var sparseResult = result as SparseMatrix;
if (sparseOther == null || sparseResult == null)
{
base.DoAdd(other, result);
return;
}
if (ReferenceEquals(this, other))
if (other is SparseMatrix sparseOther && result is SparseMatrix sparseResult)
{
if (!ReferenceEquals(this, result))
if (ReferenceEquals(this, other))
{
CopyTo(result);
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
LinearAlgebraControl.Provider.ScaleArray(2.0f, sparseResult._storage.Values, sparseResult._storage.Values);
return;
}
LinearAlgebraControl.Provider.ScaleArray(2.0f, sparseResult._storage.Values, sparseResult._storage.Values);
return;
}
SparseMatrix left;
SparseMatrix left;
if (ReferenceEquals(sparseOther, sparseResult))
{
left = this;
}
else if (ReferenceEquals(this, sparseResult))
{
left = sparseOther;
}
else
{
CopyTo(sparseResult);
left = sparseOther;
}
if (ReferenceEquals(sparseOther, sparseResult))
{
left = this;
}
else if (ReferenceEquals(this, sparseResult))
{
left = sparseOther;
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
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);
result.At(i, columnIndex, resVal);
}
}
}
else
{
CopyTo(sparseResult);
left = sparseOther;
}
var leftStorage = left._storage;
for (var i = 0; i < leftStorage.RowCount; i++)
{
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);
result.At(i, columnIndex, resVal);
}
base.DoAdd(other, result);
}
}
@ -765,59 +764,58 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <exception cref="ArgumentOutOfRangeException">If the two matrices don't have the same dimensions.</exception>
protected override void DoSubtract(Matrix<float> other, Matrix<float> result)
{
var sparseOther = other as SparseMatrix;
var sparseResult = result as SparseMatrix;
if (sparseOther == null || sparseResult == null)
{
base.DoSubtract(other, result);
return;
}
if (ReferenceEquals(this, other))
if (other is SparseMatrix sparseOther && result is SparseMatrix sparseResult)
{
result.Clear();
return;
}
if (ReferenceEquals(this, other))
{
result.Clear();
return;
}
var otherStorage = sparseOther._storage;
var otherStorage = sparseOther._storage;
if (ReferenceEquals(this, sparseResult))
{
for (var i = 0; i < otherStorage.RowCount; i++)
if (ReferenceEquals(this, sparseResult))
{
var endIndex = otherStorage.RowPointers[i + 1];
for (var j = otherStorage.RowPointers[i]; j < endIndex; j++)
for (var i = 0; i < otherStorage.RowCount; i++)
{
var columnIndex = otherStorage.ColumnIndices[j];
var resVal = sparseResult.At(i, columnIndex) - otherStorage.Values[j];
result.At(i, columnIndex, resVal);
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];
result.At(i, columnIndex, resVal);
}
}
}
}
else
{
if (!ReferenceEquals(sparseOther, sparseResult))
else
{
sparseOther.CopyTo(sparseResult);
}
if (!ReferenceEquals(sparseOther, sparseResult))
{
sparseOther.CopyTo(sparseResult);
}
sparseResult.Negate(sparseResult);
sparseResult.Negate(sparseResult);
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
var rowPointers = _storage.RowPointers;
var columnIndices = _storage.ColumnIndices;
var values = _storage.Values;
for (var i = 0; i < RowCount; i++)
{
var endIndex = rowPointers[i + 1];
for (var j = rowPointers[i]; j < endIndex; j++)
for (var i = 0; i < RowCount; i++)
{
var columnIndex = columnIndices[j];
var resVal = sparseResult.At(i, columnIndex) + values[j];
result.At(i, columnIndex, resVal);
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];
result.At(i, columnIndex, resVal);
}
}
}
}
else
{
base.DoSubtract(other, result);
}
}
/// <summary>
@ -839,8 +837,16 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return;
}
var sparseResult = result as SparseMatrix;
if (sparseResult == null)
if (result is SparseMatrix sparseResult)
{
if (!ReferenceEquals(this, result))
{
CopyTo(sparseResult);
}
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
@ -865,15 +871,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
}
}
}
else
{
if (!ReferenceEquals(this, result))
{
CopyTo(sparseResult);
}
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
}
/// <summary>
@ -1089,57 +1086,55 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">The result of the multiplication.</param>
protected override void DoTransposeAndMultiply(Matrix<float> other, Matrix<float> result)
{
var otherSparse = other as SparseMatrix;
var resultSparse = result as SparseMatrix;
if (otherSparse == null || resultSparse == null)
if (other is SparseMatrix otherSparse && result is SparseMatrix resultSparse)
{
base.DoTransposeAndMultiply(other, result);
return;
}
resultSparse.Clear();
resultSparse.Clear();
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
var otherStorage = otherSparse._storage;
for (var j = 0; j < RowCount; j++)
{
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = otherStorage.RowPointers[j + 1];
var rowPointers = _storage.RowPointers;
var values = _storage.Values;
if (startIndexOther == endIndexOther)
{
continue;
}
var otherStorage = otherSparse._storage;
for (var i = 0; i < RowCount; i++)
for (var j = 0; j < RowCount; j++)
{
// Multiply row of matrix A on row of matrix B
var startIndexThis = rowPointers[i];
var endIndexThis = rowPointers[i + 1];
var startIndexOther = otherStorage.RowPointers[j];
var endIndexOther = otherStorage.RowPointers[j + 1];
if (startIndexThis == endIndexThis)
if (startIndexOther == endIndexOther)
{
continue;
}
var sum = 0f;
for (var index = startIndexOther; index < endIndexOther; index++)
for (var i = 0; i < RowCount; i++)
{
var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
if (ind >= 0)
// Multiply row of matrix A on row of matrix B
var startIndexThis = rowPointers[i];
var endIndexThis = rowPointers[i + 1];
if (startIndexThis == endIndexThis)
{
sum += otherStorage.Values[index]*values[ind];
continue;
}
var sum = 0f;
for (var index = startIndexOther; index < endIndexOther; index++)
{
var ind = _storage.FindItem(i, otherStorage.ColumnIndices[index]);
if (ind >= 0)
{
sum += otherStorage.Values[index] * values[ind];
}
}
}
resultSparse._storage.At(i, j, sum + result.At(i, j));
resultSparse._storage.At(i, j, sum + result.At(i, j));
}
}
}
else
{
base.DoTransposeAndMultiply(other, result);
}
}
/// <summary>
@ -1266,22 +1261,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">Matrix to store the results in.</param>
protected override void DoModulus(float divisor, Matrix<float> result)
{
var sparseResult = result as SparseMatrix;
if (sparseResult == null)
if (result is SparseMatrix sparseResult)
{
base.DoModulus(divisor, result);
return;
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
var resultStorage = sparseResult._storage;
for (var index = 0; index < resultStorage.Values.Length; index++)
{
resultStorage.Values[index] = Euclid.Modulus(resultStorage.Values[index], divisor);
}
}
var resultStorage = sparseResult._storage;
for (var index = 0; index < resultStorage.Values.Length; index++)
else
{
resultStorage.Values[index] = Euclid.Modulus(resultStorage.Values[index], divisor);
base.DoModulus(divisor, result);
}
}
@ -1293,22 +1288,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">Matrix to store the results in.</param>
protected override void DoRemainder(float divisor, Matrix<float> result)
{
var sparseResult = result as SparseMatrix;
if (sparseResult == null)
if (result is SparseMatrix sparseResult)
{
base.DoRemainder(divisor, result);
return;
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
}
if (!ReferenceEquals(this, result))
{
CopyTo(result);
var resultStorage = sparseResult._storage;
for (var index = 0; index < resultStorage.Values.Length; index++)
{
resultStorage.Values[index] %= divisor;
}
}
var resultStorage = sparseResult._storage;
for (var index = 0; index < resultStorage.Values.Length; index++)
else
{
resultStorage.Values[index] %= divisor;
base.DoRemainder(divisor, result);
}
}

253
src/Numerics/LinearAlgebra/Single/SparseVector.cs

@ -192,76 +192,71 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </param>
protected override void DoAdd(Vector<float> other, Vector<float> result)
{
var otherSparse = other as SparseVector;
if (otherSparse == null)
if (other is SparseVector otherSparse && result is SparseVector resultSparse)
{
base.DoAdd(other, result);
return;
}
var resultSparse = result as SparseVector;
if (resultSparse == null)
{
base.DoAdd(other, result);
return;
}
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
{
var otherValue = otherStorage.Values[j];
if (otherValue != 0.0f)
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
{
var otherValue = otherStorage.Values[j];
if (otherValue != 0.0f)
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], otherValue);
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], otherValue);
// TODO: result can be zero, remove?
_storage.Values[i++] += otherStorage.Values[j++];
}
else
{
i++;
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] += otherStorage.Values[j++];
}
else
{
i++;
}
}
}
else
{
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
else
{
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
{
var next = _storage.Indices[i];
if (next != last)
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
{
last = next;
result.At(next, _storage.Values[i] + otherSparse.At(next));
var next = _storage.Indices[i];
if (next != last)
{
last = next;
result.At(next, _storage.Values[i] + otherSparse.At(next));
}
i++;
}
i++;
}
else
{
var next = otherStorage.Indices[j];
if (next != last)
else
{
last = next;
result.At(next, At(next) + otherStorage.Values[j]);
var next = otherStorage.Indices[j];
if (next != last)
{
last = next;
result.At(next, At(next) + otherStorage.Values[j]);
}
j++;
}
j++;
}
}
}
else
{
base.DoAdd(other, result);
}
}
/// <summary>
@ -295,76 +290,71 @@ namespace MathNet.Numerics.LinearAlgebra.Single
return;
}
var otherSparse = other as SparseVector;
if (otherSparse == null)
if (other is SparseVector otherSparse && result is SparseVector resultSparse)
{
base.DoSubtract(other, result);
return;
}
var resultSparse = result as SparseVector;
if (resultSparse == null)
{
base.DoSubtract(other, result);
return;
}
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
// TODO (ruegg, 2011-10-11): Options to optimize?
var otherStorage = otherSparse._storage;
if (ReferenceEquals(this, resultSparse))
{
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
int i = 0, j = 0;
while (j < otherStorage.ValueCount)
{
var otherValue = otherStorage.Values[j];
if (otherValue != 0.0f)
if (i >= _storage.ValueCount || _storage.Indices[i] > otherStorage.Indices[j])
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], -otherValue);
var otherValue = otherStorage.Values[j];
if (otherValue != 0.0f)
{
_storage.InsertAtIndexUnchecked(i++, otherStorage.Indices[j], -otherValue);
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] -= otherStorage.Values[j++];
}
else
{
i++;
}
j++;
}
else if (_storage.Indices[i] == otherStorage.Indices[j])
{
// TODO: result can be zero, remove?
_storage.Values[i++] -= otherStorage.Values[j++];
}
else
{
i++;
}
}
}
else
{
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
else
{
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
result.Clear();
int i = 0, j = 0, last = -1;
while (i < _storage.ValueCount || j < otherStorage.ValueCount)
{
var next = _storage.Indices[i];
if (next != last)
if (j >= otherStorage.ValueCount || i < _storage.ValueCount && _storage.Indices[i] <= otherStorage.Indices[j])
{
last = next;
result.At(next, _storage.Values[i] - otherSparse.At(next));
var next = _storage.Indices[i];
if (next != last)
{
last = next;
result.At(next, _storage.Values[i] - otherSparse.At(next));
}
i++;
}
i++;
}
else
{
var next = otherStorage.Indices[j];
if (next != last)
else
{
last = next;
result.At(next, At(next) - otherStorage.Values[j]);
var next = otherStorage.Indices[j];
if (next != last)
{
last = next;
result.At(next, At(next) - otherStorage.Values[j]);
}
j++;
}
j++;
}
}
}
else
{
base.DoSubtract(other, result);
}
}
/// <summary>
@ -373,27 +363,27 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// <param name="result">Target vector</param>
protected override void DoNegate(Vector<float> result)
{
var sparseResult = result as SparseVector;
if (sparseResult == null)
if (result is SparseVector sparseResult)
{
if (!ReferenceEquals(this, result))
{
sparseResult._storage.ValueCount = _storage.ValueCount;
sparseResult._storage.Indices = new int[_storage.ValueCount];
Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt);
sparseResult._storage.Values = new float[_storage.ValueCount];
Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
}
LinearAlgebraControl.Provider.ScaleArray(-1.0f, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], -_storage.Values[index]);
}
return;
}
if (!ReferenceEquals(this, result))
{
sparseResult._storage.ValueCount = _storage.ValueCount;
sparseResult._storage.Indices = new int[_storage.ValueCount];
Buffer.BlockCopy(_storage.Indices, 0, sparseResult._storage.Indices, 0, _storage.ValueCount * Constants.SizeOfInt);
sparseResult._storage.Values = new float[_storage.ValueCount];
Array.Copy(_storage.Values, 0, sparseResult._storage.Values, 0, _storage.ValueCount);
}
LinearAlgebraControl.Provider.ScaleArray(-1.0f, sparseResult._storage.Values, sparseResult._storage.Values);
}
/// <summary>
@ -407,16 +397,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single
/// </param>
protected override void DoMultiply(float scalar, Vector<float> result)
{
var sparseResult = result as SparseVector;
if (sparseResult == null)
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], scalar * _storage.Values[index]);
}
}
else
if (result is SparseVector sparseResult)
{
if (!ReferenceEquals(this, result))
{
@ -429,6 +410,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single
LinearAlgebraControl.Provider.ScaleArray(scalar, sparseResult._storage.Values, sparseResult._storage.Values);
}
else
{
result.Clear();
for (var index = 0; index < _storage.ValueCount; index++)
{
result.At(_storage.Indices[index], scalar * _storage.Values[index]);
}
}
}
/// <summary>

49
src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs

@ -198,43 +198,44 @@ namespace MathNet.Numerics.LinearAlgebra.Storage
return true;
}
var otherSparse = other as SparseVectorStorage<T>;
if (otherSparse == null)
if (other is SparseVectorStorage<T> otherSparse)
{
return base.Equals(other);
}
int i = 0, j = 0;
while (i < ValueCount || j < otherSparse.ValueCount)
{
if (j >= otherSparse.ValueCount || i < ValueCount && Indices[i] < otherSparse.Indices[j])
int i = 0, j = 0;
while (i < ValueCount || j < otherSparse.ValueCount)
{
if (!Zero.Equals(Values[i++]))
if (j >= otherSparse.ValueCount || i < ValueCount && Indices[i] < otherSparse.Indices[j])
{
return false;
if (!Zero.Equals(Values[i++]))
{
return false;
}
continue;
}
continue;
}
if (i >= ValueCount || j < otherSparse.ValueCount && otherSparse.Indices[j] < Indices[i])
{
if (!Zero.Equals(otherSparse.Values[j++]))
if (i >= ValueCount || j < otherSparse.ValueCount && otherSparse.Indices[j] < Indices[i])
{
if (!Zero.Equals(otherSparse.Values[j++]))
{
return false;
}
continue;
}
if (!Values[i].Equals(otherSparse.Values[j]))
{
return false;
}
continue;
}
if (!Values[i].Equals(otherSparse.Values[j]))
{
return false;
i++;
j++;
}
i++;
j++;
return true;
}
return true;
return base.Equals(other);
}
/// <summary>

Loading…
Cancel
Save