diff --git a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs index 6105bee5..1cf0f3aa 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseMatrix.cs @@ -455,21 +455,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the addition. protected override void DoAdd(Complex scalar, Matrix 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); + } } /// @@ -510,21 +510,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the subtraction. protected override void DoSubtract(Complex scalar, Matrix 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); + } } /// @@ -563,14 +563,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the multiplication. protected override void DoMultiply(Complex scalar, Matrix 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 /// The result of the multiplication. protected override void DoMultiply(Vector rightSide, Vector 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); + } } /// @@ -901,14 +897,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The matrix to store the result of the division. protected override void DoDivide(Complex divisor, Matrix 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 /// The matrix to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Matrix other, Matrix 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 /// The matrix to store the result of the pointwise division. protected override void DoPointwiseDivide(Matrix divisor, Matrix 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 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Matrix exponent, Matrix 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); } } diff --git a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs index 335fcc57..85c4078d 100644 --- a/src/Numerics/LinearAlgebra/Complex/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex/DenseVector.cs @@ -206,20 +206,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The vector to store the result of the addition. protected override void DoAdd(Complex scalar, Vector 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 /// The vector to store the result of the addition. protected override void DoAdd(Vector other, Vector 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 /// The vector to store the result of the subtraction. protected override void DoSubtract(Complex scalar, Vector 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 /// The vector to store the result of the subtraction. protected override void DoSubtract(Vector other, Vector 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 /// Target vector protected override void DoNegate(Vector 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); } /// @@ -361,14 +353,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Target vector protected override void DoConjugate(Vector 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); } /// @@ -379,14 +371,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoMultiply(Complex scalar, Vector 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); } /// @@ -396,10 +388,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The sum of a[i]*b[i] for all i. protected override Complex DoDotProduct(Vector 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); } /// @@ -409,16 +400,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The sum of conj(a[i])*b[i] for all i. protected override Complex DoConjugateDotProduct(Vector 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); } /// @@ -607,16 +601,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The vector to store the result of the pointwise division. protected override void DoPointwiseMultiply(Vector other, Vector 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 /// protected override void DoPointwiseDivide(Vector divisor, Vector 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 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Vector exponent, Vector 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); } } diff --git a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs index 2fff31e8..08319abb 100644 --- a/src/Numerics/LinearAlgebra/Complex/DiagonalMatrix.cs +++ b/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]. public override void SetDiagonal(Vector 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); } /// Calculates the induced L1 norm of this matrix. diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs index 3dec8edb..03de8c52 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseCholesky.cs @@ -94,24 +94,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -131,24 +126,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -171,19 +161,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(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."); + } } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs index d1f57bdf..7a232f96 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseGramSchmidt.cs +++ b/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); } /// @@ -182,19 +177,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs index d71a0e1d..73e33d68 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs @@ -113,24 +113,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -162,24 +157,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(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); } /// diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs index 35a9e95d..688617db 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs +++ b/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); } /// @@ -153,19 +148,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs index 352f6895..03bb4cc1 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseSvd.cs +++ b/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); } /// @@ -145,19 +140,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs index 1224f077..26fe3074 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs @@ -703,51 +703,50 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// If the two matrices don't have the same dimensions. protected override void DoAdd(Matrix other, Matrix 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); - } } /// @@ -1087,57 +1085,55 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// The result of the multiplication. protected override void DoTransposeAndMultiply(Matrix other, Matrix 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); + } } /// diff --git a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs index 5bd7e35f..6ba1fb7a 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseVector.cs @@ -191,76 +191,71 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoAdd(Vector other, Vector 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); + } } /// @@ -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); + } } /// @@ -372,27 +362,27 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// Target vector protected override void DoNegate(Vector 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); } /// @@ -434,16 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoMultiply(Complex scalar, Vector 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]); + } + } } /// diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs index 08c89197..27dab0e7 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseMatrix.cs @@ -455,21 +455,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the addition. protected override void DoAdd(Complex32 scalar, Matrix 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); + } } /// @@ -510,21 +510,21 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the subtraction. protected override void DoSubtract(Complex32 scalar, Matrix 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); + } } /// @@ -563,14 +563,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the multiplication. protected override void DoMultiply(Complex32 scalar, Matrix 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 /// The result of the multiplication. protected override void DoMultiply(Vector rightSide, Vector 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); + } } /// @@ -751,14 +747,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The result of the multiplication. protected override void DoTransposeThisAndMultiply(Vector rightSide, Vector 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); + } } /// @@ -905,14 +898,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The matrix to store the result of the division. protected override void DoDivide(Complex32 divisor, Matrix 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 /// The matrix to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Matrix other, Matrix 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 /// The matrix to store the result of the pointwise division. protected override void DoPointwiseDivide(Matrix divisor, Matrix 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 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Matrix exponent, Matrix 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); } } diff --git a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs index 5644df8b..a244003f 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/DenseVector.cs @@ -206,20 +206,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The vector to store the result of the addition. protected override void DoAdd(Complex32 scalar, Vector 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 /// The vector to store the result of the addition. protected override void DoAdd(Vector other, Vector 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 /// The vector to store the result of the subtraction. protected override void DoSubtract(Complex32 scalar, Vector 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 /// The vector to store the result of the subtraction. protected override void DoSubtract(Vector other, Vector 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 /// Target vector protected override void DoNegate(Vector 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); } /// @@ -361,14 +353,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// Target vector protected override void DoConjugate(Vector 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); } /// @@ -379,14 +371,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoMultiply(Complex32 scalar, Vector 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); } /// @@ -396,10 +388,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The sum of a[i]*b[i] for all i. protected override Complex32 DoDotProduct(Vector 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); } /// @@ -409,16 +400,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The sum of conj(a[i])*b[i] for all i. protected override Complex32 DoConjugateDotProduct(Vector 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); } /// @@ -607,16 +601,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The vector to store the result of the pointwise division. protected override void DoPointwiseMultiply(Vector other, Vector 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 /// protected override void DoPointwiseDivide(Vector divisor, Vector 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 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Vector exponent, Vector 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); } } diff --git a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs index 6f866966..b8dee3d2 100644 --- a/src/Numerics/LinearAlgebra/Complex32/DiagonalMatrix.cs +++ b/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]. public override void SetDiagonal(Vector 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); } /// Calculates the induced L1 norm of this matrix. diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs index 77973b03..6b84e096 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseCholesky.cs @@ -94,24 +94,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -131,24 +126,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -171,19 +161,20 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(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."); + } } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs index 71b1ef98..afd641cf 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseGramSchmidt.cs +++ b/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); } /// @@ -182,19 +177,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs index aa1d6a19..856ed21f 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs @@ -113,24 +113,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -162,24 +157,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(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); } /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs index 28355089..4921136a 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs +++ b/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); } /// @@ -153,19 +148,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs index 9aa73b30..cb00ce16 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseSvd.cs +++ b/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); } /// @@ -145,19 +140,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs index 3902f255..dab55e11 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs @@ -703,51 +703,50 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// If the two matrices don't have the same dimensions. protected override void DoAdd(Matrix other, Matrix 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 /// If the two matrices don't have the same dimensions. protected override void DoSubtract(Matrix other, Matrix 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); + } } /// @@ -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); - } } /// @@ -1086,57 +1083,55 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// The result of the multiplication. protected override void DoTransposeAndMultiply(Matrix other, Matrix 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); + } } /// diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs index bcf89326..f25c8e96 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseVector.cs @@ -191,76 +191,71 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoAdd(Vector other, Vector 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); + } } /// @@ -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); + } } /// @@ -372,27 +362,27 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// Target vector protected override void DoNegate(Vector 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); } /// @@ -434,16 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoMultiply(Complex32 scalar, Vector 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]); + } + } } /// diff --git a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs index bbf65c08..24665233 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs @@ -438,21 +438,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the addition. protected override void DoAdd(double scalar, Matrix 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); + } } /// @@ -493,21 +493,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the subtraction. protected override void DoSubtract(double scalar, Matrix 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); + } } /// @@ -546,14 +546,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the multiplication. protected override void DoMultiply(double scalar, Matrix 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 /// The result of the multiplication. protected override void DoMultiply(Vector rightSide, Vector 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); + } } /// @@ -681,14 +677,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The result of the multiplication. protected override void DoTransposeThisAndMultiply(Vector rightSide, Vector 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); + } } /// @@ -760,14 +753,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The matrix to store the result of the division. protected override void DoDivide(double divisor, Matrix 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 /// The matrix to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Matrix other, Matrix 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 /// The matrix to store the result of the pointwise division. protected override void DoPointwiseDivide(Matrix divisor, Matrix 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 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Matrix exponent, Matrix 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 /// Matrix to store the results in. protected override void DoModulus(double divisor, Matrix 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); + } } /// @@ -869,21 +852,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoModulusByThis(double dividend, Matrix 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); + } } /// @@ -894,26 +877,26 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Matrix to store the results in. protected override void DoRemainder(double divisor, Matrix 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); + } } /// @@ -924,21 +907,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoRemainderByThis(double dividend, Matrix 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); + } } /// diff --git a/src/Numerics/LinearAlgebra/Double/DenseVector.cs b/src/Numerics/LinearAlgebra/Double/DenseVector.cs index fa1cec37..80ab7346 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs @@ -206,20 +206,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The vector to store the result of the addition. protected override void DoAdd(double scalar, Vector 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 /// The vector to store the result of the addition. protected override void DoAdd(Vector other, Vector 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 /// The vector to store the result of the subtraction. protected override void DoSubtract(double scalar, Vector 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 /// The vector to store the result of the subtraction. protected override void DoSubtract(Vector other, Vector 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 /// Target vector protected override void DoNegate(Vector 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); } /// @@ -373,14 +365,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoMultiply(double scalar, Vector 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); } /// @@ -390,10 +382,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The sum of a[i]*b[i] for all i. protected override double DoDotProduct(Vector 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); } /// @@ -473,12 +464,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoModulus(double divisor, Vector 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); + } } /// @@ -498,21 +488,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoRemainder(double divisor, Vector 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); + } } /// @@ -689,16 +678,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The vector to store the result of the pointwise division. protected override void DoPointwiseMultiply(Vector other, Vector 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 /// protected override void DoPointwiseDivide(Vector divisor, Vector 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 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Vector exponent, Vector 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); } } diff --git a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs index eea86d64..3af9b702 100644 --- a/src/Numerics/LinearAlgebra/Double/DiagonalMatrix.cs +++ b/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]. public override void SetDiagonal(Vector 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); } /// Calculates the induced L1 norm of this matrix. @@ -843,21 +842,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Matrix to store the results in. protected override void DoModulus(double divisor, Matrix 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); + } } /// @@ -868,21 +867,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoModulusByThis(double dividend, Matrix 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); + } } /// @@ -893,21 +892,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Matrix to store the results in. protected override void DoRemainder(double divisor, Matrix 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); + } } /// @@ -918,21 +917,21 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// A vector to store the results in. protected override void DoRemainderByThis(double dividend, Matrix 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); + } } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs index 795b6498..d3a75e27 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseCholesky.cs @@ -92,24 +92,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -129,24 +124,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -169,19 +159,20 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(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."); + } } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs index fcd9fd1b..7adba9d9 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseGramSchmidt.cs +++ b/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); } /// @@ -180,19 +175,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs index b6aaef81..9cbade23 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs @@ -111,24 +111,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -160,24 +155,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(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); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs index 079ed421..fec6f37a 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs +++ b/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); } /// @@ -151,19 +146,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs index b667f6a0..13539043 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs +++ b/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); } /// @@ -143,19 +138,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs index f2b3d1cd..7b48ec96 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs @@ -703,51 +703,50 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// If the two matrices don't have the same dimensions. protected override void DoAdd(Matrix other, Matrix 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 /// If the two matrices don't have the same dimensions. protected override void DoSubtract(Matrix other, Matrix 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); + } } /// @@ -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); - } } /// @@ -1086,55 +1082,53 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// The result of the multiplication. protected override void DoTransposeAndMultiply(Matrix other, Matrix 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); + } } /// @@ -1261,22 +1255,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Matrix to store the results in. protected override void DoModulus(double divisor, Matrix 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 /// Matrix to store the results in. protected override void DoRemainder(double divisor, Matrix 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); } } diff --git a/src/Numerics/LinearAlgebra/Double/SparseVector.cs b/src/Numerics/LinearAlgebra/Double/SparseVector.cs index 7cf97629..1c710405 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseVector.cs @@ -191,76 +191,71 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoAdd(Vector other, Vector 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); + } } /// @@ -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); + } } /// @@ -372,27 +362,27 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// Target vector protected override void DoNegate(Vector 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); } /// @@ -406,16 +396,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoMultiply(double scalar, Vector 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]); + } + } } /// diff --git a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs index 5e5f5964..f17c1ed8 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseMatrix.cs @@ -438,21 +438,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the addition. protected override void DoAdd(float scalar, Matrix 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); + } } /// @@ -493,21 +493,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the subtraction. protected override void DoSubtract(float scalar, Matrix 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); + } } /// @@ -546,14 +546,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the multiplication. protected override void DoMultiply(float scalar, Matrix 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 /// The result of the multiplication. protected override void DoMultiply(Vector rightSide, Vector 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); + } } /// @@ -757,14 +753,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The matrix to store the result of the division. protected override void DoDivide(float divisor, Matrix 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 /// The matrix to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Matrix other, Matrix 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 /// The matrix to store the result of the pointwise division. protected override void DoPointwiseDivide(Matrix divisor, Matrix 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 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Matrix exponent, Matrix 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 /// Matrix to store the results in. protected override void DoModulus(float divisor, Matrix 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); + } } /// @@ -866,21 +852,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoModulusByThis(float dividend, Matrix 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); + } } /// @@ -891,26 +877,26 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Matrix to store the results in. protected override void DoRemainder(float divisor, Matrix 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); + } } /// @@ -921,21 +907,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoRemainderByThis(float dividend, Matrix 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); + } } /// diff --git a/src/Numerics/LinearAlgebra/Single/DenseVector.cs b/src/Numerics/LinearAlgebra/Single/DenseVector.cs index a8a789f6..8bbb5cce 100644 --- a/src/Numerics/LinearAlgebra/Single/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Single/DenseVector.cs @@ -206,20 +206,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The vector to store the result of the addition. protected override void DoAdd(float scalar, Vector 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 /// The vector to store the result of the addition. protected override void DoAdd(Vector other, Vector 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 /// The vector to store the result of the subtraction. protected override void DoSubtract(float scalar, Vector 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 /// The vector to store the result of the subtraction. protected override void DoSubtract(Vector other, Vector 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 /// Target vector protected override void DoNegate(Vector 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); } /// @@ -363,14 +355,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoMultiply(float scalar, Vector 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); } /// @@ -380,10 +372,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The sum of a[i]*b[i] for all i. protected override float DoDotProduct(Vector 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); } /// @@ -463,12 +454,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoModulus(float divisor, Vector 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); + } } /// @@ -488,21 +478,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoRemainder(float divisor, Vector 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); + } } /// @@ -680,16 +669,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The vector to store the result of the pointwise multiplication. protected override void DoPointwiseMultiply(Vector other, Vector 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 /// protected override void DoPointwiseDivide(Vector divisor, Vector 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 /// The vector to store the result of the pointwise power. protected override void DoPointwisePower(Vector exponent, Vector 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); } } diff --git a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs b/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs index e16dce4a..0deaabd1 100644 --- a/src/Numerics/LinearAlgebra/Single/DiagonalMatrix.cs +++ b/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]. public override void SetDiagonal(Vector 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); } /// Calculates the induced L1 norm of this matrix. @@ -843,21 +842,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Matrix to store the results in. protected override void DoModulus(float divisor, Matrix 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); + } } /// @@ -868,21 +867,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoModulusByThis(float dividend, Matrix 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); + } } /// @@ -893,21 +892,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Matrix to store the results in. protected override void DoRemainder(float divisor, Matrix 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); + } } /// @@ -918,21 +917,21 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// A vector to store the results in. protected override void DoRemainderByThis(float dividend, Matrix 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); + } } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs index f86e68a9..5ead5410 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseCholesky.cs @@ -92,24 +92,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -129,24 +124,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -169,19 +159,20 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(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."); + } } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs index 853f8fa1..b6fe7a7c 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseGramSchmidt.cs +++ b/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); } /// @@ -180,19 +175,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs index fb407307..c3e24945 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs @@ -111,24 +111,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(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); } /// @@ -160,24 +155,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(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); } /// diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs index c0f7785e..8250e9b1 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs +++ b/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); } /// @@ -151,19 +146,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs index f8f9d51a..24ca5177 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseSvd.cs +++ b/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); } /// @@ -144,19 +139,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization throw Matrix.DimensionsDontMatch(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); } } } diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs index e0f92453..bd44c4f2 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs @@ -708,51 +708,50 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// If the two matrices don't have the same dimensions. protected override void DoAdd(Matrix other, Matrix 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 /// If the two matrices don't have the same dimensions. protected override void DoSubtract(Matrix other, Matrix 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); + } } /// @@ -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); - } } /// @@ -1089,57 +1086,55 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// The result of the multiplication. protected override void DoTransposeAndMultiply(Matrix other, Matrix 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); + } } /// @@ -1266,22 +1261,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Matrix to store the results in. protected override void DoModulus(float divisor, Matrix 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 /// Matrix to store the results in. protected override void DoRemainder(float divisor, Matrix 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); } } diff --git a/src/Numerics/LinearAlgebra/Single/SparseVector.cs b/src/Numerics/LinearAlgebra/Single/SparseVector.cs index 6105a1e0..0ef59b64 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseVector.cs @@ -192,76 +192,71 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoAdd(Vector other, Vector 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); + } } /// @@ -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); + } } /// @@ -373,27 +363,27 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// Target vector protected override void DoNegate(Vector 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); } /// @@ -407,16 +397,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoMultiply(float scalar, Vector 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]); + } + } } /// diff --git a/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs b/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs index f8c19a3e..51fb8f23 100644 --- a/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs +++ b/src/Numerics/LinearAlgebra/Storage/SparseVectorStorage.cs @@ -198,43 +198,44 @@ namespace MathNet.Numerics.LinearAlgebra.Storage return true; } - var otherSparse = other as SparseVectorStorage; - if (otherSparse == null) + if (other is SparseVectorStorage 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); } ///