diff --git a/src/Numerics/LinearAlgebra/Double/DenseVector.cs b/src/Numerics/LinearAlgebra/Double/DenseVector.cs index d45772a3..84ca992f 100644 --- a/src/Numerics/LinearAlgebra/Double/DenseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/DenseVector.cs @@ -29,6 +29,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double using System; using System.Collections.Generic; using System.Globalization; + using System.Linq; using Distributions; using NumberTheory; using Properties; @@ -303,10 +304,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 0.0) { - return this.Clone(); + return Clone(); } - var copy = (DenseVector) this.Clone(); + var copy = (DenseVector)Clone(); CommonParallel.For( 0, Data.Length, @@ -343,7 +344,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double CommonParallel.For( 0, Data.Length, - index => dense.Data[index] = this.Data[index] + scalar); + index => dense.Data[index] = Data[index] + scalar); } } @@ -372,12 +373,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { return base.Add(other); } - else - { - var copy = (DenseVector)this.Clone(); - Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, 1.0, denseVector.Data); - return copy; - } + + var copy = (DenseVector)Clone(); + Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, 1.0, denseVector.Data); + return copy; } /// @@ -408,7 +407,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - var tmp = this.Add(other); + var tmp = Add(other); tmp.CopyTo(result); } else @@ -425,7 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double CommonParallel.For( 0, Data.Length, - index => result[index] = this.Data[index] + other[index]); + index => result[index] = Data[index] + other[index]); } } } @@ -484,10 +483,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 0.0) { - return this.Clone(); + return Clone(); } - var copy = (DenseVector)this.Clone(); + var copy = (DenseVector)Clone(); CommonParallel.For( 0, Data.Length, @@ -517,14 +516,14 @@ namespace MathNet.Numerics.LinearAlgebra.Double var dense = result as DenseVector; if (dense == null) { - base.Add(scalar, result); + base.Subtract(scalar, result); } else { CommonParallel.For( 0, Data.Length, - index => dense.Data[index] = this.Data[index] - scalar); + index => dense.Data[index] = Data[index] - scalar); } } @@ -553,12 +552,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { return base.Subtract(other); } - else - { - var copy = (DenseVector)this.Clone(); - Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, -1.0, denseVector.Data); - return copy; - } + + var copy = (DenseVector)Clone(); + Control.LinearAlgebraProvider.AddVectorToScaledVector(copy.Data, -1.0, denseVector.Data); + return copy; } /// @@ -589,8 +586,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - var tmp = result.CreateVector(result.Count); - Subtract(other, tmp); + var tmp = Subtract(other); tmp.CopyTo(result); } else @@ -607,7 +603,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double CommonParallel.For( 0, Data.Length, - index => result[index] = this.Data[index] - other[index]); + index => result[index] = Data[index] - other[index]); } } } @@ -681,10 +677,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 1.0) { - return this.Clone(); + return Clone(); } - var copy = (DenseVector)this.Clone(); + var copy = (DenseVector)Clone(); Control.LinearAlgebraProvider.ScaleArray(scalar, copy.Data); return copy; } @@ -732,7 +728,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentNullException("leftSide"); } - return (DenseVector) leftSide.Multiply(rightSide); + return (DenseVector)leftSide.Multiply(rightSide); } /// @@ -749,7 +745,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentNullException("rightSide"); } - return (DenseVector) rightSide.Multiply(leftSide); + return (DenseVector)rightSide.Multiply(leftSide); } /// @@ -794,7 +790,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentNullException("leftSide"); } - return (DenseVector) leftSide.Multiply(1.0 / rightSide); + return (DenseVector)leftSide.Multiply(1.0 / rightSide); } /// @@ -1013,15 +1009,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double { return base.PointwiseMultiply(other); } - else - { - var copy = (DenseVector)this.Clone(); - CommonParallel.For( - 0, - Count, - index => copy[index] *= other[index]); - return copy; - } + + var copy = (DenseVector)Clone(); + CommonParallel.For( + 0, + Count, + index => copy[index] *= other[index]); + return copy; } /// @@ -1057,8 +1051,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - var tmp = result.CreateVector(result.Count); - PointwiseMultiply(other, tmp); + var tmp = PointwiseMultiply(other); tmp.CopyTo(result); } else @@ -1073,7 +1066,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double CommonParallel.For( 0, Data.Length, - index => dense.Data[index] = this.Data[index] * other[index]); + index => dense.Data[index] = Data[index] * other[index]); } } } @@ -1103,15 +1096,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double { return base.PointwiseMultiply(other); } - else - { - var copy = (DenseVector)this.Clone(); - CommonParallel.For( - 0, - Count, - index => copy[index] /= other[index]); - return copy; - } + + var copy = (DenseVector)Clone(); + CommonParallel.For( + 0, + Count, + index => copy[index] /= other[index]); + return copy; } /// @@ -1147,7 +1138,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - var tmp = this.PointwiseDivide(other); + var tmp = PointwiseDivide(other); tmp.CopyTo(result); } else @@ -1162,7 +1153,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double CommonParallel.For( 0, Data.Length, - index => dense.Data[index] = this.Data[index] / other[index]); + index => dense.Data[index] = Data[index] / other[index]); } } } @@ -1279,41 +1270,35 @@ namespace MathNet.Numerics.LinearAlgebra.Double { throw new ArgumentOutOfRangeException("p"); } - else if (1.0 == p) + + if (1.0 == p) { return CommonParallel.Aggregate( 0, Count, index => Math.Abs(Data[index])); } - else if (2.0 == p) - { - var sum = 0.0; - for (var i = 0; i < Data.Length; i++) - { - sum = SpecialFunctions.Hypotenuse(sum, Data[i]); - } - - return sum; + if (2.0 == p) + { + return Data.Aggregate(0.0, SpecialFunctions.Hypotenuse); } - else if (Double.IsPositiveInfinity(p)) + + if (Double.IsPositiveInfinity(p)) { return CommonParallel.Select( 0, Count, - (index, localData) => localData = Math.Max(localData, Math.Abs(Data[index])), + (index, localData) => Math.Max(localData, Math.Abs(Data[index])), Math.Max); } - else - { - var sum = CommonParallel.Aggregate( - 0, - Count, - index => Math.Pow(Math.Abs(Data[index]), p)); - return Math.Pow(sum, 1.0 / p); - } + var sum = CommonParallel.Aggregate( + 0, + Count, + index => Math.Pow(Math.Abs(Data[index]), p)); + + return Math.Pow(sum, 1.0 / p); } #endregion diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs index db81cfac..a2ea9a3e 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseSvd.cs @@ -165,13 +165,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization var dinput = input as DenseVector; if (dinput == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); } var dresult = result as DenseVector; if (dresult == null) { - throw new NotImplementedException("Can only do QR factorization for dense vectors at the moment."); + throw new NotImplementedException("Can only do SVD factorization for dense vectors at the moment."); } Control.LinearAlgebraProvider.SvdSolveFactored(MatrixU.RowCount, MatrixVT.ColumnCount, ((DenseVector)VectorS).Data, ((DenseMatrix)MatrixU).Data, ((DenseMatrix)MatrixVT).Data, dinput.Data, 1, dresult.Data); diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs index c13ac636..1b73a71c 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/QR.cs @@ -97,7 +97,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization } /// - /// Gets the determinant of the matrix for which the QR matrix was computed. + /// Gets the absolute determinant value of the matrix for which the QR matrix was computed. /// public virtual double Determinant { diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs index c3279cbc..76c5c3bf 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs @@ -290,8 +290,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative Vector vecSdash = new DenseVector(residuals.Count); Vector t = new DenseVector(residuals.Count); - Vector mult = new DenseVector(residuals.Count); - // create some temporary double variables that are needed // to hold values in between iterations double currentRho = 0; @@ -319,11 +317,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative var beta = (currentRho / oldRho) * (alpha / omega); // p_i = r_(i-1) + beta_(i-1)(p_(i-1) - omega_(i-1) * nu_(i-1)) - nu.Multiply(-omega, mult); - vecP.Add(mult); - - vecP.Multiply(beta); - vecP.Add(residuals); + vecP.Add(nu.Multiply(-omega), vecP); + + vecP.Multiply(beta, vecP); + vecP.Add(residuals, vecP); } else { @@ -341,8 +338,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative alpha = currentRho * 1 / tempResiduals.DotProduct(nu); // s = r_(i-1) - alpha_i nu_i - nu.Multiply(-alpha, mult); - residuals.Add(mult, vecS); + residuals.Add(nu.Multiply(-alpha), vecS); // Check if we're converged. If so then stop. Otherwise continue; // Calculate the temporary result. @@ -350,8 +346,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative // temp. Others will be used in the calculation later on. // x_i = x_(i-1) + alpha_i * p^_i + s^_i vecPdash.Multiply(alpha, temp); - temp.Add(vecSdash); - temp.Add(result); + temp.Add(vecSdash, temp); + temp.Add(result, temp); // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, temp, input, vecS)) @@ -383,14 +379,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative omega = t.DotProduct(vecS) / t.DotProduct(t); // x_i = x_(i-1) + alpha_i p^ + omega_i s^ - vecSdash.Multiply(omega, mult); - result.Add(mult); + result.Add(vecSdash.Multiply(omega), result); - vecPdash.Multiply(alpha, mult); - result.Add(mult); + //vecPdash.Multiply(alpha, vecPdash); + result.Add(vecPdash.Multiply(alpha), result); t.Multiply(-omega, residuals); - residuals.Add(vecS); + residuals.Add(vecS, residuals); // for continuation it is necessary that omega_i != 0.0 // If omega is only 1 ULP from zero then we fail. @@ -427,10 +422,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative matrix.Multiply(x, residual); // Do not use residual = residual.Negate() because it creates another object - residual.Multiply(-1); + residual.Multiply(-1, residual); // residual + b - residual.Add(b); + residual.Add(b, residual); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs index d981caf6..b783170a 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs @@ -359,7 +359,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative Vector z = new DenseVector(residuals.Count); Vector temp = new DenseVector(residuals.Count); - Vector mult = new DenseVector(residuals.Count); // for (k = 0, 1, .... ) var iterationNumber = 0; @@ -368,8 +367,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative // p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1)) p.Subtract(u, temp); - temp.Multiply(beta, mult); - residuals.Add(mult, p); + residuals.Add(temp.Multiply(beta), p); // Solve M b_k = p_k _preconditioner.Approximate(p, temp); @@ -384,15 +382,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative s.Subtract(w, temp); t.Subtract(residuals, y); - temp.Multiply(alpha, mult); - y.Add(mult); + y.Add(temp.Multiply(alpha), y); // Store the old value of t in t0 t.CopyTo(t0); // t_k = r_k - alpha_k s_k - s.Multiply(-alpha, mult); - residuals.Add(mult, t); + residuals.Add(s.Multiply(-alpha), t); // Solve M d_k = t_k _preconditioner.Approximate(t, temp); @@ -450,48 +446,39 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative } // u_k = sigma_k s_k + eta_k (t_(k-1) - r_k + beta_(k-1) u_(k-1)) - u.Multiply(beta, mult); - t0.Add(mult, temp); - - temp.Subtract(residuals); - temp.Multiply(eta); + t0.Add(u.Multiply(beta), temp); - s.Multiply(sigma, mult); - temp.Add(mult, u); + temp.Subtract(residuals, temp); + temp.Multiply(eta, temp); - // z_k = sigma_k r_k +_ eta_k z_(k-1) - alpha_k u_k - z.Multiply(eta); + temp.Add(s.Multiply(sigma), u); - u.Multiply(-alpha, mult); - z.Add(mult); + // z_k = sigma_k r_k +_ eta_k z_(k-1) - alpha_k u_k + z.Multiply(eta, z); + z.Add(u.Multiply(-alpha), z); - residuals.Multiply(sigma, mult); - z.Add(mult); + z.Add(residuals.Multiply(sigma), z); // x_(k+1) = x_k + alpha_k p_k + z_k - p.Multiply(alpha, mult); - xtemp.Add(mult); - - xtemp.Add(z); + xtemp.Add(p.Multiply(alpha), xtemp); + + xtemp.Add(z, xtemp); // r_(k+1) = t_k - eta_k y_k - sigma_k c_k // Copy the old residuals to a temp vector because we'll // need those in the next step residuals.CopyTo(t0); - y.Multiply(-eta, mult); - t.Add(mult, residuals); + t.Add(y.Multiply(-eta), residuals); - c.Multiply(-sigma, mult); - residuals.Add(mult); + residuals.Add(c.Multiply(-sigma), residuals); // beta_k = alpha_k / sigma_k * (r*_0 * r_(k+1)) / (r*_0 * r_k) // But first we check if there is a possible NaN. If so just reset beta to zero. beta = (!sigma.AlmostEqual(0, 1)) ? alpha / sigma * rdash.DotProduct(residuals) / rdash.DotProduct(t0) : 0; // w_k = c_k + beta_k s_k - s.Multiply(beta, mult); - c.Add(mult, w); + c.Add(s.Multiply(beta), w); // Get the real value _preconditioner.Approximate(xtemp, result); @@ -520,10 +507,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative { // -Ax = residual matrix.Multiply(x, residual); - residual.Multiply(-1); + residual.Multiply(-1, residual); // residual + b - residual.Add(b); + residual.Add(b, residual); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs index cd9f5d29..0f8ccadf 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs @@ -395,8 +395,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative Vector zg = new DenseVector(residuals.Count); Vector zw = new DenseVector(residuals.Count); - Vector mult = new DenseVector(residuals.Count); - var d = CreateVectorArray(_startingVectors.Count, residuals.Count); // g_0 = r_0 @@ -426,8 +424,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative var alpha = _startingVectors[0].DotProduct(residuals) / c[k - 1]; // u_(jk+1) = r_((j-1)k+k) - alpha_(jk+1) w_((j-1)k+k) - w[k - 1].Multiply(-alpha, mult); - residuals.Add(mult, u); + residuals.Add(w[k - 1].Multiply(-alpha), u); // SOLVE M u~_(jk+1) = u_(jk+1) _preconditioner.Approximate(u, temp1); @@ -448,18 +445,17 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative rho = -u.DotProduct(temp) / rho; // x_(jk+1) = x_((j-1)k_k) - rho_(j+1) u~_(jk+1) + alpha_(jk+1) g~_((j-1)k+k) - utemp.Multiply(-rho, mult); - xtemp.Add(mult); + xtemp.Add(utemp.Multiply(-rho), xtemp); - gtemp.Multiply(alpha); - xtemp.Add(gtemp); + gtemp.Multiply(alpha, gtemp); + xtemp.Add(gtemp, xtemp); // r_(jk+1) = rho_(j+1) A u~_(jk+1) + u_(jk+1) u.CopyTo(residuals); // Reuse temp - temp.Multiply(rho); - residuals.Add(temp); + temp.Multiply(rho, temp); + residuals.Add(temp, residuals); // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) @@ -498,16 +494,13 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative beta = -_startingVectors[s + 1].DotProduct(zd) / c[s]; // z_d = z_d + beta^(jk+i)_((j-1)k+s) d_((j-1)k+s) - d[s].Multiply(beta, mult); - zd.Add(mult); + zd.Add(d[s].Multiply(beta), zd); // z_g = z_g + beta^(jk+i)_((j-1)k+s) g_((j-1)k+s) - g[s].Multiply(beta, mult); - zg.Add(mult); + zg.Add(g[s].Multiply(beta), zg); // z_w = z_w + beta^(jk+i)_((j-1)k+s) w_((j-1)k+s) - w[s].Multiply(beta, mult); - zw.Add(mult); + zw.Add(w[s].Multiply(beta), zw); } } @@ -518,18 +511,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative } // beta^(jk+i)_((j-1)k+k) = -(q^T_1 (r_(jk+1) + rho_(j+1) z_w)) / (rho_(j+1) c_((j-1)k+k)) - zw.Multiply(rho, mult); - residuals.Add(mult, temp); + residuals.Add(zw.Multiply(rho), temp); beta = -_startingVectors[0].DotProduct(temp) / beta; // z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k) - g[k - 1].Multiply(beta, mult); - zg.Add(mult); + zg.Add(g[k - 1].Multiply(beta), zg); // z_w = rho_(j+1) (z_w + beta^(jk+i)_((j-1)k+k) w_((j-1)k+k)) - w[k - 1].Multiply(beta, mult); - zw.Add(mult); - zw.Multiply(rho); + zw.Add(w[k - 1].Multiply(beta), zw); + zw.Multiply(rho, zw); // z_d = r_(jk+i) + z_w residuals.Add(zw, zd); @@ -541,12 +531,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative beta = -_startingVectors[s + 1].DotProduct(zd) / c[s]; // z_d = z_d + beta^(jk+i)_(jk+s) * d_(jk+s) - d[s].Multiply(beta, mult); - zd.Add(mult); + zd.Add(d[s].Multiply(beta), zd); // z_g = z_g + beta^(jk+i)_(jk+s) * g_(jk+s) - g[s].Multiply(beta, mult); - zg.Add(mult); + zg.Add(g[s].Multiply(beta), zg); } // d_(jk+i) = z_d - u_(jk+i) @@ -569,22 +557,19 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative alpha = _startingVectors[i + 1].DotProduct(u) / c[i]; // u_(jk+i+1) = u_(jk+i) - alpha_(jk+i+1) d_(jk+i) - d[i].Multiply(-alpha, mult); - u.Add(mult); + u.Add(d[i].Multiply(-alpha), u); // SOLVE M g~_(jk+i) = g_(jk+i) _preconditioner.Approximate(g[i], gtemp); // x_(jk+i+1) = x_(jk+i) + rho_(j+1) alpha_(jk+i+1) g~_(jk+i) - gtemp.Multiply(rho * alpha, mult); - xtemp.Add(mult); + xtemp.Add(gtemp.Multiply(rho * alpha), xtemp); // w_(jk+i) = A g~_(jk+i) matrix.Multiply(gtemp, w[i]); // r_(jk+i+1) = r_(jk+i) - rho_(j+1) alpha_(jk+i+1) w_(jk+i) - w[i].Multiply(-rho * alpha, mult); - residuals.Add(mult); + residuals.Add(w[i].Multiply(-rho * alpha), residuals); // We can check the residuals here if they're close if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) @@ -657,7 +642,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative result.Add(orthogonalMatrix.Column(i)); // Normalize the result vector - result[i].Multiply(1 / result[i].Norm(2)); + result[i].Multiply(1 / result[i].Norm(2), result[i]); } return result; @@ -691,10 +676,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative { // -Ax = residual matrix.Multiply(x, residual); - residual.Multiply(-1); + residual.Multiply(-1, residual); // residual + b - residual.Add(b); + residual.Add(b, residual); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs index b06df5a9..78d01fa0 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs @@ -273,8 +273,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative // Temp vectors var temp = new DenseVector(input.Count); - var mult = new DenseVector(input.Count); - + // Initialize var startNorm = input.Norm(2); @@ -316,8 +315,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative alpha = rho / sigma; // yOdd = yEven - alpha * v - v.Multiply(-alpha, mult); - yeven.Add(mult, yodd); + yeven.Add(v.Multiply(-alpha), yodd); // Solve M temp = yOdd _preconditioner.Approximate(yodd, temp); @@ -333,8 +331,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative var yinternal = IsEven(iterationNumber) ? yeven : yodd; // pseudoResiduals = pseudoResiduals - alpha * uOdd - uinternal.Multiply(-alpha, mult); - pseudoResiduals.Add(mult); + pseudoResiduals.Add(uinternal.Multiply(-alpha), pseudoResiduals); // d = yOdd + theta * theta * eta / alpha * d d.Multiply(theta * theta * eta / alpha, temp); @@ -351,8 +348,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative eta = c * c * alpha; // x = x + eta * d - d.Multiply(eta, mult); - x.Add(mult); + x.Add(d.Multiply(eta), x); // Check convergence and see if we can bail if (!ShouldContinue(iterationNumber, result, input, pseudoResiduals)) @@ -390,8 +386,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative rho = rhoNew; // yOdd = pseudoResiduals + beta * yOdd - yodd.Multiply(beta, mult); - pseudoResiduals.Add(mult, yeven); + pseudoResiduals.Add(yodd.Multiply(beta), yeven); // Solve M temp = yOdd _preconditioner.Approximate(yeven, temp); @@ -400,11 +395,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative matrix.Multiply(temp, ueven); // v = uEven + beta * (uOdd + beta * v) - v.Multiply(beta, mult); - uodd.Add(mult, temp); + uodd.Add(v.Multiply(beta), temp); - temp.Multiply(beta, mult); - ueven.Add(mult, v); + ueven.Add(temp.Multiply(beta), v); } // Calculate the real values @@ -425,10 +418,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative { // -Ax = residual matrix.Multiply(x, residual); - residual.Multiply(-1); + residual.Multiply(-1, residual); // residual + b - residual.Add(b); + residual.Add(b, residual); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/Ilutp.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/Ilutp.cs index 8938e19b..9c8d65b2 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/Ilutp.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Preconditioners/Ilutp.cs @@ -426,8 +426,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Preconditioners rowVector[k] = 0.0; } - rowVector.Multiply(workVector[j]); - workVector.Subtract(rowVector); + rowVector.Multiply(workVector[j], rowVector); + workVector.Subtract(rowVector, workVector); } } } diff --git a/src/Numerics/LinearAlgebra/Double/SparseVector.cs b/src/Numerics/LinearAlgebra/Double/SparseVector.cs index b1a24145..238cbe09 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseVector.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseVector.cs @@ -351,14 +351,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 0.0) { - return this.Clone(); + return Clone(); } - var copy = (SparseVector)this.Clone(); + var copy = (SparseVector)Clone(); for (var i = 0; i < Count; i++) { copy[i] += scalar; } + return copy; } @@ -409,12 +410,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { return base.Add(other); } - else - { - var copy = (SparseVector)this.Clone(); - copy.AddScaledSparseVector(1.0, sparseVector); - return copy; - } + + var copy = (SparseVector)Clone(); + copy.AddScaledSparseVector(1.0, sparseVector); + return copy; } /// @@ -514,7 +513,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - var tmp = this.Add(other); + var tmp = Add(other); tmp.CopyTo(result); } else @@ -593,14 +592,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 0.0) { - return this.Clone(); + return Clone(); } - var copy = (SparseVector)this.Clone(); + var copy = (SparseVector)Clone(); for (var i = 0; i < Count; i++) { copy[i] -= scalar; } + return copy; } @@ -651,12 +651,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { return base.Subtract(other); } - else - { - var copy = (SparseVector)this.Clone(); - copy.AddScaledSparseVector(-1.0, sparseVector); - return copy; - } + + var copy = (SparseVector)Clone(); + copy.AddScaledSparseVector(-1.0, sparseVector); + return copy; } /// @@ -687,7 +685,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - var tmp = this.Subtract(other); + var tmp = Subtract(other); tmp.CopyTo(result); } else @@ -789,18 +787,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 1.0) { - return this.Clone(); + return Clone(); } - else if (scalar == 0) + + if (scalar == 0) { - var copy = this.Clone(); + var copy = Clone(); copy.Clear(); // Set array empty return copy; } else { - - var copy = (SparseVector)this.Clone(); + var copy = (SparseVector)Clone(); Control.LinearAlgebraProvider.ScaleArray(scalar, copy._nonZeroValues); return copy; } @@ -912,7 +910,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentNullException("leftSide"); } - return (SparseVector) leftSide.Multiply(1.0 / rightSide); + return (SparseVector)leftSide.Multiply(1.0 / rightSide); } /// @@ -1103,12 +1101,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double } var copy = new SparseVector(Count); - for (int i = 0; i < this._nonZeroIndices.Length; i++) + for (int i = 0; i < _nonZeroIndices.Length; i++) { - var d = this._nonZeroValues[i] * other[this._nonZeroIndices[i]]; + var d = _nonZeroValues[i] * other[_nonZeroIndices[i]]; if (d != 0.0) { - copy[this._nonZeroIndices[i]] = d; + copy[_nonZeroIndices[i]] = d; } } @@ -1238,7 +1236,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (Double.IsPositiveInfinity(p)) { - return CommonParallel.Select(0, NonZerosCount, (index, localData) => localData = Math.Max(localData, Math.Abs(_nonZeroValues[index])), Math.Max); + return CommonParallel.Select(0, NonZerosCount, (index, localData) => Math.Max(localData, Math.Abs(_nonZeroValues[index])), Math.Max); } var sum = CommonParallel.Aggregate( diff --git a/src/Numerics/LinearAlgebra/Double/Vector.cs b/src/Numerics/LinearAlgebra/Double/Vector.cs index 6d6085e5..38064bab 100644 --- a/src/Numerics/LinearAlgebra/Double/Vector.cs +++ b/src/Numerics/LinearAlgebra/Double/Vector.cs @@ -125,10 +125,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 0.0) { - return this.Clone(); + return Clone(); } - var copy = this.Clone(); + var copy = Clone(); CommonParallel.For( 0, Count, @@ -185,7 +185,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// public virtual Vector Plus() { - return this.Clone(); + return Clone(); } /// @@ -213,7 +213,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); } - var copy = this.Clone(); + var copy = Clone(); CommonParallel.For( 0, Count, @@ -256,7 +256,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - var tmp = this.Add(other); + var tmp = Add(other); tmp.CopyTo(result); } else @@ -279,10 +279,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 0.0) { - return this.Clone(); + return Clone(); } - var copy = this.Clone(); + var copy = Clone(); CommonParallel.For( 0, Count, @@ -367,7 +367,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); } - var copy = this.Clone(); + var copy = Clone(); CommonParallel.For( 0, Count, @@ -434,10 +434,10 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 1.0) { - return this.Clone(); + return Clone(); } - var copy = this.Clone(); + var copy = Clone(); CommonParallel.For( 0, Count, @@ -530,7 +530,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double { if (scalar == 1.0) { - return this.Clone(); + return Clone(); } return Multiply(1.0 / scalar); @@ -593,7 +593,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); } - var copy = this.Clone(); + var copy = Clone(); CommonParallel.For( 0, Count, @@ -634,7 +634,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - var tmp = this.PointwiseMultiply(other); + var tmp = PointwiseMultiply(other); tmp.CopyTo(result); } else @@ -665,7 +665,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); } - var copy = this.Clone(); + var copy = Clone(); CommonParallel.For( 0, Count, @@ -706,7 +706,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double if (ReferenceEquals(this, result) || ReferenceEquals(other, result)) { - var tmp = this.PointwiseDivide(other); + var tmp = PointwiseDivide(other); tmp.CopyTo(result); } else @@ -1156,23 +1156,22 @@ namespace MathNet.Numerics.LinearAlgebra.Double { throw new ArgumentOutOfRangeException("p"); } - else if (Double.IsPositiveInfinity(p)) + + if (Double.IsPositiveInfinity(p)) { return CommonParallel.Select( 0, Count, - (index, localData) => localData = Math.Max(localData, Math.Abs(this[index])), + (index, localData) => Math.Max(localData, Math.Abs(this[index])), Math.Max); } - else - { - var sum = CommonParallel.Aggregate( - 0, - Count, - index => Math.Pow(Math.Abs(this[index]), p)); - return Math.Pow(sum, 1.0 / p); - } + var sum = CommonParallel.Aggregate( + 0, + Count, + index => Math.Pow(Math.Abs(this[index]), p)); + + return Math.Pow(sum, 1.0 / p); } /// @@ -1192,7 +1191,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double } var norm = Norm(p); - var clone = this.Clone(); + var clone = Clone(); if (norm == 0.0) { return clone; diff --git a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs index a405cc3c..6714e40d 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/MatrixTests.Arithmetic.cs @@ -695,7 +695,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Row(-4, ExpectedException = typeof(ArgumentOutOfRangeException))] public void NormalizeColumns(int pValue) { - var matrix = this.TestMatrices["Singular3x3"]; + var matrix = this.TestMatrices["Square4x4"]; var result = matrix.NormalizeColumns(pValue); for (var j = 0; j < result.ColumnCount; j++) { @@ -710,7 +710,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double [Row(-3, ExpectedException = typeof(ArgumentOutOfRangeException))] public void NormalizeRows(int pValue) { - var matrix = this.TestMatrices["Singular3x3"].NormalizeRows(pValue); + var matrix = this.TestMatrices["Square4x4"].NormalizeRows(pValue); for (var i = 0; i < matrix.RowCount; i++) { var row = matrix.Row(i); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs index 1e851a84..c980fda5 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs @@ -40,11 +40,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative public void SolveUnitMatrixAndBackMultiply() { // Create the identity matrix - Matrix matrix = new SparseMatrix(100); - for (var i = 0; i < matrix.RowCount; i++) - { - matrix[i, i] = 1.0; - } + Matrix matrix = SparseMatrix.Identity(100); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); @@ -84,11 +80,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative public void SolveScaledUnitMatrixAndBackMultiply() { // Create the identity matrix - Matrix matrix = new SparseMatrix(100); - for (var i = 0; i < matrix.RowCount; i++) - { - matrix[i, i] = 1.0; - } + Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number matrix.Multiply(System.Math.PI); @@ -201,5 +193,69 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative Assert.IsTrue(System.Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#05-" + i); } } + + [Test] + [Row(4)] + [Row(8)] + [Row(10)] + [MultipleAsserts] + public void CanSolveForRandomVector(int order) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); + var vectorb = MatrixLoader.GenerateRandomDenseVector(order); + + var monitor = new Iterator(new IIterationStopCriterium[] + { + new IterationCountStopCriterium(1000), + new ResidualStopCriterium(1e-10), + }); + var solver = new BiCgStab(monitor); + + var resultx = solver.Solve(matrixA, vectorb); + Assert.AreEqual(matrixA.ColumnCount, resultx.Count); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < order; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-7); + } + } + + [Test] + [Row(4)] + [Row(8)] + [Row(10)] + [MultipleAsserts] + public void CanSolveForRandomMatrix(int order) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); + + var monitor = new Iterator(new IIterationStopCriterium[] + { + new IterationCountStopCriterium(1000), + new ResidualStopCriterium(1e-10) + }); + var solver = new BiCgStab(monitor); + var matrixX = solver.Solve(matrixA, matrixB); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-7); + } + } + } } } diff --git a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs index a4154c7b..a56236fc 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs @@ -40,11 +40,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative public void SolveUnitMatrixAndBackMultiply() { // Create the identity matrix - Matrix matrix = new SparseMatrix(100); - for (var i = 0; i < matrix.RowCount; i++) - { - matrix[i, i] = 1.0; - } + Matrix matrix = SparseMatrix.Identity(100); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); @@ -84,11 +80,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative public void SolveScaledUnitMatrixAndBackMultiply() { // Create the identity matrix - Matrix matrix = new SparseMatrix(100); - for (var i = 0; i < matrix.RowCount; i++) - { - matrix[i, i] = 1.0; - } + Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number matrix.Multiply(System.Math.PI); @@ -203,5 +195,69 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative Assert.IsTrue(System.Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#05-" + i); } } + + [Test] + [Row(4)] + [Row(8)] + [Row(10)] + [MultipleAsserts] + public void CanSolveForRandomVector(int order) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); + var vectorb = MatrixLoader.GenerateRandomDenseVector(order); + + var monitor = new Iterator(new IIterationStopCriterium[] + { + new IterationCountStopCriterium(1000), + new ResidualStopCriterium(1e-10), + }); + var solver = new GpBiCg(monitor); + + var resultx = solver.Solve(matrixA, vectorb); + Assert.AreEqual(matrixA.ColumnCount, resultx.Count); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < order; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-7); + } + } + + [Test] + [Row(4)] + [Row(8)] + [Row(10)] + [MultipleAsserts] + public void CanSolveForRandomMatrix(int order) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); + + var monitor = new Iterator(new IIterationStopCriterium[] + { + new IterationCountStopCriterium(1000), + new ResidualStopCriterium(1e-10) + }); + var solver = new GpBiCg(monitor); + var matrixX = solver.Solve(matrixA, matrixB); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-7); + } + } + } } } \ No newline at end of file diff --git a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs index 92883217..3582e146 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs @@ -40,11 +40,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative public void SolveUnitMatrixAndBackMultiply() { // Create the identity matrix - Matrix matrix = new SparseMatrix(100); - for (var i = 0; i < matrix.RowCount; i++) - { - matrix[i, i] = 1.0; - } + Matrix matrix = SparseMatrix.Identity(100); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); @@ -85,11 +81,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative public void SolveScaledUnitMatrixAndBackMultiply() { // Create the identity matrix - Matrix matrix = new SparseMatrix(100); - for (var i = 0; i < matrix.RowCount; i++) - { - matrix[i, i] = 1.0; - } + Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number matrix.Multiply(System.Math.PI); @@ -201,5 +193,69 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative Assert.IsTrue(System.Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#05-" + i); } } + + [Test] + [Row(4)] + [Row(8)] + [Row(10)] + [MultipleAsserts] + public void CanSolveForRandomVector(int order) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); + var vectorb = MatrixLoader.GenerateRandomDenseVector(order); + + var monitor = new Iterator(new IIterationStopCriterium[] + { + new IterationCountStopCriterium(1000), + new ResidualStopCriterium(1e-10), + }); + var solver = new MlkBiCgStab(monitor); + + var resultx = solver.Solve(matrixA, vectorb); + Assert.AreEqual(matrixA.ColumnCount, resultx.Count); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < order; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-7); + } + } + + [Test] + [Row(4)] + [Row(8)] + [Row(10)] + [MultipleAsserts] + public void CanSolveForRandomMatrix(int order) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); + + var monitor = new Iterator(new IIterationStopCriterium[] + { + new IterationCountStopCriterium(1000), + new ResidualStopCriterium(1e-10) + }); + var solver = new MlkBiCgStab(monitor); + var matrixX = solver.Solve(matrixA, matrixB); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-7); + } + } + } } } diff --git a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs index 20267552..0b123eb7 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs @@ -40,11 +40,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative public void SolveUnitMatrixAndBackMultiply() { // Create the identity matrix - Matrix matrix = new SparseMatrix(100); - for (var i = 0; i < matrix.RowCount; i++) - { - matrix[i, i] = 1.0; - } + Matrix matrix = SparseMatrix.Identity(100); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); @@ -85,11 +81,8 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative public void SolveScaledUnitMatrixAndBackMultiply() { // Create the identity matrix - Matrix matrix = new SparseMatrix(100); - for (var i = 0; i < matrix.RowCount; i++) - { - matrix[i, i] = 1.0; - } + Matrix matrix = SparseMatrix.Identity(100); + // Scale it with a funny number matrix.Multiply(System.Math.PI); @@ -200,5 +193,69 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative Assert.IsTrue(System.Math.Abs(y[i] - z[i]).IsSmaller(ConvergenceBoundary, 1), "#05-" + i); } } + + [Test] + [Row(4)] + [Row(8)] + [Row(10)] + [MultipleAsserts] + public void CanSolveForRandomVector(int order) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); + var vectorb = MatrixLoader.GenerateRandomDenseVector(order); + + var monitor = new Iterator(new IIterationStopCriterium[] + { + new IterationCountStopCriterium(1000), + new ResidualStopCriterium(1e-10), + }); + var solver = new TFQMR(monitor); + + var resultx = solver.Solve(matrixA, vectorb); + Assert.AreEqual(matrixA.ColumnCount, resultx.Count); + + var bReconstruct = matrixA * resultx; + + // Check the reconstruction. + for (var i = 0; i < order; i++) + { + Assert.AreApproximatelyEqual(vectorb[i], bReconstruct[i], 1e-7); + } + } + + [Test] + [Row(4)] + [Row(8)] + [Row(10)] + [MultipleAsserts] + public void CanSolveForRandomMatrix(int order) + { + var matrixA = MatrixLoader.GenerateRandomDenseMatrix(order, order); + var matrixB = MatrixLoader.GenerateRandomDenseMatrix(order, order); + + var monitor = new Iterator(new IIterationStopCriterium[] + { + new IterationCountStopCriterium(1000), + new ResidualStopCriterium(1e-10) + }); + var solver = new TFQMR(monitor); + var matrixX = solver.Solve(matrixA, matrixB); + + // The solution X row dimension is equal to the column dimension of A + Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount); + // The solution X has the same number of columns as B + Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount); + + var matrixBReconstruct = matrixA * matrixX; + + // Check the reconstruction. + for (var i = 0; i < matrixB.RowCount; i++) + { + for (var j = 0; j < matrixB.ColumnCount; j++) + { + Assert.AreApproximatelyEqual(matrixB[i, j], matrixBReconstruct[i, j], 1.0e-7); + } + } + } } }