From 9bf22ea0111e4ace5a9b56ebf41f2cbb41b01b0f Mon Sep 17 00:00:00 2001 From: Marcus Cuda Date: Sun, 21 Nov 2010 19:49:01 +0800 Subject: [PATCH] opt: merged Andriy's changes to reuse vector and matrices instead of creating new ones --- .../Complex/Solvers/Iterative/BiCgStab.cs | 47 +++++++++---- .../Complex/Solvers/Iterative/GpBiCg.cs | 59 +++++++++++----- .../Complex/Solvers/Iterative/MlkBiCgStab.cs | 70 ++++++++++++++----- .../Complex/Solvers/Iterative/TFQMR.cs | 30 ++++++-- .../LinearAlgebra/Complex/SparseMatrix.cs | 25 +++++++ src/Numerics/LinearAlgebra/Complex/Vector.cs | 1 - .../Complex32/Solvers/Iterative/BiCgStab.cs | 47 +++++++++---- .../Complex32/Solvers/Iterative/GpBiCg.cs | 59 +++++++++++----- .../Solvers/Iterative/MlkBiCgStab.cs | 70 ++++++++++++++----- .../Complex32/Solvers/Iterative/TFQMR.cs | 30 ++++++-- .../LinearAlgebra/Complex32/SparseMatrix.cs | 25 +++++++ .../LinearAlgebra/Complex32/Vector.cs | 1 - src/Numerics/LinearAlgebra/Double/Matrix.cs | 2 +- .../Double/Solvers/Iterative/BiCgStab.cs | 49 ++++++++----- .../Double/Solvers/Iterative/GpBiCg.cs | 62 +++++++++++----- .../Double/Solvers/Iterative/MlkBiCgStab.cs | 70 ++++++++++++++----- .../Double/Solvers/Iterative/TFQMR.cs | 27 +++++-- .../LinearAlgebra/Double/SparseMatrix.cs | 25 +++++++ src/Numerics/LinearAlgebra/Double/Vector.cs | 1 - .../Single/Solvers/Iterative/BiCgStab.cs | 47 +++++++++---- .../Single/Solvers/Iterative/GpBiCg.cs | 58 ++++++++++----- .../Single/Solvers/Iterative/MlkBiCgStab.cs | 70 ++++++++++++++----- .../Single/Solvers/Iterative/TFQMR.cs | 29 ++++++-- .../LinearAlgebra/Single/SparseMatrix.cs | 25 +++++++ src/Numerics/LinearAlgebra/Single/Vector.cs | 3 - .../Complex/Solvers/Iterative/BiCgStabTest.cs | 2 +- .../Complex/Solvers/Iterative/GpBiCgTest.cs | 2 +- .../Solvers/Iterative/MlkBiCgStabTest.cs | 2 +- .../Complex/Solvers/Iterative/TFQMRTest.cs | 2 +- .../Solvers/Iterative/BiCgStabTest.cs | 2 +- .../Complex32/Solvers/Iterative/GpBiCgTest.cs | 2 +- .../Solvers/Iterative/MlkBiCgStabTest.cs | 2 +- .../Complex32/Solvers/Iterative/TFQMRTest.cs | 2 +- .../Double/Solvers/Iterative/BiCgStabTest.cs | 2 +- .../Double/Solvers/Iterative/GpBiCgTest.cs | 2 +- .../Solvers/Iterative/MlkBiCgStabTest.cs | 2 +- .../Double/Solvers/Iterative/TFQMRTest.cs | 2 +- .../Single/Solvers/Iterative/BiCgStabTest.cs | 2 +- .../Single/Solvers/Iterative/GpBiCgTest.cs | 2 +- .../Solvers/Iterative/MlkBiCgStabTest.cs | 2 +- .../Single/Solvers/Iterative/TFQMRTest.cs | 2 +- 41 files changed, 707 insertions(+), 257 deletions(-) diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs index 356eefca..ca2f2d06 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs @@ -259,7 +259,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - + + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -282,9 +287,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative // Choose r~ (for example, r~ = r_0) var tempResiduals = residuals.Clone(); - var temp = result.Clone(); - // create five temporary vectors needed to hold temporary + // create seven temporary vectors needed to hold temporary // coefficients. All vectors are mangled in each iteration. // These are defined here to prevent stressing the garbage collector Vector vecP = new DenseVector(residuals.Count); @@ -292,7 +296,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative Vector nu = new DenseVector(residuals.Count); Vector vecS = new DenseVector(residuals.Count); Vector vecSdash = new DenseVector(residuals.Count); - Vector t = new DenseVector(residuals.Count); + Vector temp = new DenseVector(residuals.Count); + Vector temp2 = new DenseVector(residuals.Count); // create some temporary double variables that are needed // to hold values in between iterations @@ -321,10 +326,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.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)) - vecP.Add(nu.Multiply(-omega), vecP); + nu.Multiply(-omega, temp); + vecP.Add(temp, temp2); + temp2.CopyTo(vecP); vecP.Multiply(beta, vecP); - vecP.Add(residuals, vecP); + vecP.Add(residuals, temp2); + temp2.CopyTo(vecP); } else { @@ -342,7 +350,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative alpha = currentRho * 1 / tempResiduals.DotProduct(nu); // s = r_(i-1) - alpha_i nu_i - residuals.Add(nu.Multiply(-alpha), vecS); + nu.Multiply(-alpha, temp); + residuals.Add(temp, vecS); // Check if we're converged. If so then stop. Otherwise continue; // Calculate the temporary result. @@ -350,8 +359,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.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); - temp.Add(result, temp); + temp.Add(vecSdash, temp2); + temp2.CopyTo(temp); + temp.Add(result, temp2); + temp2.CopyTo(temp); // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, temp, input, vecS)) @@ -377,17 +388,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative _preconditioner.Approximate(vecS, vecSdash); // temp = As~ - matrix.Multiply(vecSdash, t); + matrix.Multiply(vecSdash, temp); // omega_i = temp^T s / temp^T temp - omega = t.DotProduct(vecS) / t.DotProduct(t); + omega = temp.DotProduct(vecS) / temp.DotProduct(temp); // x_i = x_(i-1) + alpha_i p^ + omega_i s^ - result.Add(vecSdash.Multiply(omega), result); - result.Add(vecPdash.Multiply(alpha), result); + temp.Multiply(-omega, residuals); + residuals.Add(vecS, temp2); + temp2.CopyTo(residuals); + + vecSdash.Multiply(omega, temp); + result.Add(temp, temp2); + temp2.CopyTo(result); - t.Multiply(-omega, residuals); - residuals.Add(vecS, residuals); + vecPdash.Multiply(alpha, temp); + result.Add(temp, temp2); + temp2.CopyTo(result); // for continuation it is necessary that omega_i != 0.0 // If omega is only 1 ULP from zero then we fail. diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs index 2b98d81e..58c70cea 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs @@ -315,8 +315,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - // Initialize the solver fields + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) { @@ -363,15 +367,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative Vector z = new DenseVector(residuals.Count); Vector temp = new DenseVector(residuals.Count); - + Vector temp2 = new DenseVector(residuals.Count); + Vector temp3 = new DenseVector(residuals.Count); + // for (k = 0, 1, .... ) var iterationNumber = 0; while (ShouldContinue(iterationNumber, xtemp, input, residuals)) { // p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1)) p.Subtract(u, temp); - - residuals.Add(temp.Multiply(beta), p); + + temp.Multiply(beta, temp2); + residuals.Add(temp2, p); // Solve M b_k = p_k _preconditioner.Approximate(p, temp); @@ -385,14 +392,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative // y_k = t_(k-1) - r_k - alpha_k * w_(k-1) + alpha_k s_k s.Subtract(w, temp); t.Subtract(residuals, y); - - y.Add(temp.Multiply(alpha), y); + + temp.Multiply(alpha, temp2); + y.Add(temp2, temp3); + temp3.CopyTo(y); // Store the old value of t in t0 t.CopyTo(t0); // t_k = r_k - alpha_k s_k - residuals.Add(s.Multiply(-alpha), t); + s.Multiply(-alpha, temp2); + residuals.Add(temp2, t); // Solve M d_k = t_k _preconditioner.Approximate(t, temp); @@ -450,38 +460,53 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative } // u_k = sigma_k s_k + eta_k (t_(k-1) - r_k + beta_(k-1) u_(k-1)) - t0.Add(u.Multiply(beta), temp); + u.Multiply(beta, temp2); + t0.Add(temp2, temp); - temp.Subtract(residuals, temp); + temp.Subtract(residuals, temp3); + temp3.CopyTo(temp); temp.Multiply(eta, temp); - temp.Add(s.Multiply(sigma), u); + s.Multiply(sigma, temp2); + temp.Add(temp2, u); // 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); - z.Add(residuals.Multiply(sigma), z); + u.Multiply(-alpha, temp2); + z.Add(temp2, temp3); + temp3.CopyTo(z); + + residuals.Multiply(sigma, temp2); + z.Add(temp2, temp3); + temp3.CopyTo(z); // x_(k+1) = x_k + alpha_k p_k + z_k - xtemp.Add(p.Multiply(alpha), xtemp); + p.Multiply(alpha, temp2); + xtemp.Add(temp2, temp3); + temp3.CopyTo(xtemp); - xtemp.Add(z, xtemp); + xtemp.Add(z, temp3); + temp3.CopyTo(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); - t.Add(y.Multiply(-eta), residuals); + y.Multiply(-eta, temp2); + t.Add(temp2, residuals); - residuals.Add(c.Multiply(-sigma), residuals); + c.Multiply(-sigma, temp2); + residuals.Add(temp2, temp3); + temp3.CopyTo(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.Real.AlmostEqual(0, 1) || !sigma.Imaginary.AlmostEqual(0, 1)) ? alpha / sigma * rdash.DotProduct(residuals) / rdash.DotProduct(t0) : 0; // w_k = c_k + beta_k s_k - c.Add(s.Multiply(beta), w); + s.Multiply(beta, temp2); + c.Add(temp2, w); // Get the real value _preconditioner.Approximate(xtemp, result); diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs index e3621bb7..7671217d 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs @@ -336,6 +336,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -394,6 +399,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative Vector utemp = new DenseVector(residuals.Count); Vector temp = new DenseVector(residuals.Count); Vector temp1 = new DenseVector(residuals.Count); + Vector temp2 = new DenseVector(residuals.Count); Vector zd = new DenseVector(residuals.Count); Vector zg = new DenseVector(residuals.Count); @@ -428,7 +434,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.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) - residuals.Add(w[k - 1].Multiply(-alpha), u); + w[k - 1].Multiply(-alpha, temp); + residuals.Add(temp, u); // SOLVE M u~_(jk+1) = u_(jk+1) _preconditioner.Approximate(u, temp1); @@ -448,18 +455,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.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) - xtemp.Add(utemp.Multiply(-rho), xtemp); - - 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, temp); - residuals.Add(temp, residuals); + residuals.Add(temp, temp2); + temp2.CopyTo(residuals); + + // 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, temp); + xtemp.Add(temp, temp2); + temp2.CopyTo(xtemp); + + gtemp.Multiply(alpha, gtemp); + xtemp.Add(gtemp, temp2); + temp2.CopyTo(xtemp); // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) @@ -498,13 +509,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.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) - zd.Add(d[s].Multiply(beta), zd); + d[s].Multiply(beta, temp); + zd.Add(temp, temp2); + temp2.CopyTo(zd); // z_g = z_g + beta^(jk+i)_((j-1)k+s) g_((j-1)k+s) - zg.Add(g[s].Multiply(beta), zg); + g[s].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); // z_w = z_w + beta^(jk+i)_((j-1)k+s) w_((j-1)k+s) - zw.Add(w[s].Multiply(beta), zw); + w[s].Multiply(beta, temp); + zw.Add(temp, temp2); + temp2.CopyTo(zw); } } @@ -515,14 +532,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.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)) - residuals.Add(zw.Multiply(rho), temp); + zw.Multiply(rho, temp2); + residuals.Add(temp2, temp); beta = -_startingVectors[0].DotProduct(temp) / beta; // z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k) - zg.Add(g[k - 1].Multiply(beta), zg); + g[k - 1].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); // z_w = rho_(j+1) (z_w + beta^(jk+i)_((j-1)k+k) w_((j-1)k+k)) - zw.Add(w[k - 1].Multiply(beta), zw); + w[k - 1].Multiply(beta, temp); + zw.Add(temp, temp2); + temp2.CopyTo(zw); zw.Multiply(rho, zw); // z_d = r_(jk+i) + z_w @@ -535,10 +557,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative beta = -_startingVectors[s + 1].DotProduct(zd) / c[s]; // z_d = z_d + beta^(jk+i)_(jk+s) * d_(jk+s) - zd.Add(d[s].Multiply(beta), zd); + d[s].Multiply(beta, temp); + zd.Add(temp, temp2); + temp2.CopyTo(zd); // z_g = z_g + beta^(jk+i)_(jk+s) * g_(jk+s) - zg.Add(g[s].Multiply(beta), zg); + g[s].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); } // d_(jk+i) = z_d - u_(jk+i) @@ -561,19 +587,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative alpha = _startingVectors[i + 1].DotProduct(u) / c[i]; // u_(jk+i+1) = u_(jk+i) - alpha_(jk+i+1) d_(jk+i) - u.Add(d[i].Multiply(-alpha), u); + d[i].Multiply(-alpha, temp); + u.Add(temp, temp2); + temp2.CopyTo(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) - xtemp.Add(gtemp.Multiply(rho * alpha), xtemp); + gtemp.Multiply(rho * alpha, temp); + xtemp.Add(temp, temp2); + temp2.CopyTo(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) - residuals.Add(w[i].Multiply(-rho * alpha), residuals); + w[i].Multiply(-rho * alpha, temp); + residuals.Add(temp, temp2); + temp2.CopyTo(residuals); // We can check the residuals here if they're close if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs index 6cd59fde..b8efc344 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs @@ -248,6 +248,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -277,7 +282,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative // Temp vectors var temp = new DenseVector(input.Count); - + var temp1 = new DenseVector(input.Count); + var temp2 = new DenseVector(input.Count); + // Initialize var startNorm = input.Norm(2); @@ -319,7 +326,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative alpha = rho / sigma; // yOdd = yEven - alpha * v - yeven.Add(v.Multiply(-alpha), yodd); + v.Multiply(-alpha, temp1); + yeven.Add(temp1, yodd); // Solve M temp = yOdd _preconditioner.Approximate(yodd, temp); @@ -335,7 +343,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative var yinternal = IsEven(iterationNumber) ? yeven : yodd; // pseudoResiduals = pseudoResiduals - alpha * uOdd - pseudoResiduals.Add(uinternal.Multiply(-alpha), pseudoResiduals); + uinternal.Multiply(-alpha, temp1); + pseudoResiduals.Add(temp1, temp2); + temp2.CopyTo(pseudoResiduals); // d = yOdd + theta * theta * eta / alpha * d d.Multiply(theta * theta * eta / alpha, temp); @@ -352,7 +362,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative eta = c * c * alpha; // x = x + eta * d - x.Add(d.Multiply(eta), x); + d.Multiply(eta, temp1); + x.Add(temp1, temp2); + temp2.CopyTo(x); // Check convergence and see if we can bail if (!ShouldContinue(iterationNumber, result, input, pseudoResiduals)) @@ -390,7 +402,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative rho = rhoNew; // yOdd = pseudoResiduals + beta * yOdd - pseudoResiduals.Add(yodd.Multiply(beta), yeven); + yodd.Multiply(beta, temp1); + pseudoResiduals.Add(temp1, yeven); // Solve M temp = yOdd _preconditioner.Approximate(yeven, temp); @@ -399,8 +412,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative matrix.Multiply(temp, ueven); // v = uEven + beta * (uOdd + beta * v) - uodd.Add(v.Multiply(beta), temp); - ueven.Add(temp.Multiply(beta), v); + v.Multiply(beta, temp1); + uodd.Add(temp1, temp); + + temp.Multiply(beta, temp1); + ueven.Add(temp1, v); } // Calculate the real values diff --git a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs index d20a9c3c..4ba63ff3 100644 --- a/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex/SparseMatrix.cs @@ -1282,6 +1282,31 @@ namespace MathNet.Numerics.LinearAlgebra.Complex } } + /// + /// Multiplies this matrix with a vector and places the results into the result vector. + /// + /// The vector to multiply with. + /// The result of the multiplication. + protected override void DoMultiply(Vector rightSide, Vector result) + { + for (var row = 0; row < RowCount; row++) + { + // Get the begin / end index for the current row + var startIndex = _rowIndex[row]; + var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + if (startIndex == endIndex) + { + continue; + } + + var sum = CommonParallel.Aggregate( + startIndex, + endIndex, + index => _nonZeroValues[index] * rightSide[_columnIndices[index]]); + result[row] = sum; + } + } + /// /// Multiplies this matrix with transpose of another matrix and places the results into the result matrix. /// diff --git a/src/Numerics/LinearAlgebra/Complex/Vector.cs b/src/Numerics/LinearAlgebra/Complex/Vector.cs index ffe9b52e..e0440e25 100644 --- a/src/Numerics/LinearAlgebra/Complex/Vector.cs +++ b/src/Numerics/LinearAlgebra/Complex/Vector.cs @@ -111,7 +111,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex /// protected override void DoSubtract(Vector other, Vector result) { - CopyTo(result); CommonParallel.For( 0, Count, diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs index e09f30ed..3f832752 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs @@ -259,7 +259,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - + + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -282,9 +287,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative // Choose r~ (for example, r~ = r_0) var tempResiduals = residuals.Clone(); - var temp = result.Clone(); - // create five temporary vectors needed to hold temporary + // create seven temporary vectors needed to hold temporary // coefficients. All vectors are mangled in each iteration. // These are defined here to prevent stressing the garbage collector Vector vecP = new DenseVector(residuals.Count); @@ -292,7 +296,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative Vector nu = new DenseVector(residuals.Count); Vector vecS = new DenseVector(residuals.Count); Vector vecSdash = new DenseVector(residuals.Count); - Vector t = new DenseVector(residuals.Count); + Vector temp = new DenseVector(residuals.Count); + Vector temp2 = new DenseVector(residuals.Count); // create some temporary float variables that are needed // to hold values in between iterations @@ -321,10 +326,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.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)) - vecP.Add(nu.Multiply(-omega), vecP); + nu.Multiply(-omega, temp); + vecP.Add(temp, temp2); + temp2.CopyTo(vecP); vecP.Multiply(beta, vecP); - vecP.Add(residuals, vecP); + vecP.Add(residuals, temp2); + temp2.CopyTo(vecP); } else { @@ -342,7 +350,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative alpha = currentRho * 1 / tempResiduals.DotProduct(nu); // s = r_(i-1) - alpha_i nu_i - residuals.Add(nu.Multiply(-alpha), vecS); + nu.Multiply(-alpha, temp); + residuals.Add(temp, vecS); // Check if we're converged. If so then stop. Otherwise continue; // Calculate the temporary result. @@ -350,8 +359,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.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); - temp.Add(result, temp); + temp.Add(vecSdash, temp2); + temp2.CopyTo(temp); + temp.Add(result, temp2); + temp2.CopyTo(temp); // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, temp, input, vecS)) @@ -377,17 +388,23 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative _preconditioner.Approximate(vecS, vecSdash); // temp = As~ - matrix.Multiply(vecSdash, t); + matrix.Multiply(vecSdash, temp); // omega_i = temp^T s / temp^T temp - omega = t.DotProduct(vecS) / t.DotProduct(t); + omega = temp.DotProduct(vecS) / temp.DotProduct(temp); // x_i = x_(i-1) + alpha_i p^ + omega_i s^ - result.Add(vecSdash.Multiply(omega), result); - result.Add(vecPdash.Multiply(alpha), result); + temp.Multiply(-omega, residuals); + residuals.Add(vecS, temp2); + temp2.CopyTo(residuals); + + vecSdash.Multiply(omega, temp); + result.Add(temp, temp2); + temp2.CopyTo(result); - t.Multiply(-omega, residuals); - residuals.Add(vecS, residuals); + vecPdash.Multiply(alpha, temp); + result.Add(temp, temp2); + temp2.CopyTo(result); // for continuation it is necessary that omega_i != 0.0f // If omega is only 1 ULP from zero then we fail. diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs index e6cac702..26e35eb8 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs @@ -315,8 +315,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - // Initialize the solver fields + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) { @@ -363,15 +367,18 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative Vector z = new DenseVector(residuals.Count); Vector temp = new DenseVector(residuals.Count); - + Vector temp2 = new DenseVector(residuals.Count); + Vector temp3 = new DenseVector(residuals.Count); + // for (k = 0, 1, .... ) var iterationNumber = 0; while (ShouldContinue(iterationNumber, xtemp, input, residuals)) { // p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1)) p.Subtract(u, temp); - - residuals.Add(temp.Multiply(beta), p); + + temp.Multiply(beta, temp2); + residuals.Add(temp2, p); // Solve M b_k = p_k _preconditioner.Approximate(p, temp); @@ -385,14 +392,17 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative // y_k = t_(k-1) - r_k - alpha_k * w_(k-1) + alpha_k s_k s.Subtract(w, temp); t.Subtract(residuals, y); - - y.Add(temp.Multiply(alpha), y); + + temp.Multiply(alpha, temp2); + y.Add(temp2, temp3); + temp3.CopyTo(y); // Store the old value of t in t0 t.CopyTo(t0); // t_k = r_k - alpha_k s_k - residuals.Add(s.Multiply(-alpha), t); + s.Multiply(-alpha, temp2); + residuals.Add(temp2, t); // Solve M d_k = t_k _preconditioner.Approximate(t, temp); @@ -450,38 +460,53 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative } // u_k = sigma_k s_k + eta_k (t_(k-1) - r_k + beta_(k-1) u_(k-1)) - t0.Add(u.Multiply(beta), temp); + u.Multiply(beta, temp2); + t0.Add(temp2, temp); - temp.Subtract(residuals, temp); + temp.Subtract(residuals, temp3); + temp3.CopyTo(temp); temp.Multiply(eta, temp); - temp.Add(s.Multiply(sigma), u); + s.Multiply(sigma, temp2); + temp.Add(temp2, u); // 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); - z.Add(residuals.Multiply(sigma), z); + u.Multiply(-alpha, temp2); + z.Add(temp2, temp3); + temp3.CopyTo(z); + + residuals.Multiply(sigma, temp2); + z.Add(temp2, temp3); + temp3.CopyTo(z); // x_(k+1) = x_k + alpha_k p_k + z_k - xtemp.Add(p.Multiply(alpha), xtemp); + p.Multiply(alpha, temp2); + xtemp.Add(temp2, temp3); + temp3.CopyTo(xtemp); - xtemp.Add(z, xtemp); + xtemp.Add(z, temp3); + temp3.CopyTo(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); - t.Add(y.Multiply(-eta), residuals); + y.Multiply(-eta, temp2); + t.Add(temp2, residuals); - residuals.Add(c.Multiply(-sigma), residuals); + c.Multiply(-sigma, temp2); + residuals.Add(temp2, temp3); + temp3.CopyTo(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.Real.AlmostEqual(0, 1) || !sigma.Imaginary.AlmostEqual(0, 1)) ? alpha / sigma * rdash.DotProduct(residuals) / rdash.DotProduct(t0) : 0; // w_k = c_k + beta_k s_k - c.Add(s.Multiply(beta), w); + s.Multiply(beta, temp2); + c.Add(temp2, w); // Get the real value _preconditioner.Approximate(xtemp, result); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs index 0b6fbf19..81ee45fb 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs @@ -336,6 +336,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -394,6 +399,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative Vector utemp = new DenseVector(residuals.Count); Vector temp = new DenseVector(residuals.Count); Vector temp1 = new DenseVector(residuals.Count); + Vector temp2 = new DenseVector(residuals.Count); Vector zd = new DenseVector(residuals.Count); Vector zg = new DenseVector(residuals.Count); @@ -428,7 +434,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.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) - residuals.Add(w[k - 1].Multiply(-alpha), u); + w[k - 1].Multiply(-alpha, temp); + residuals.Add(temp, u); // SOLVE M u~_(jk+1) = u_(jk+1) _preconditioner.Approximate(u, temp1); @@ -448,18 +455,22 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.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) - xtemp.Add(utemp.Multiply(-rho), xtemp); - - 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, temp); - residuals.Add(temp, residuals); + residuals.Add(temp, temp2); + temp2.CopyTo(residuals); + + // 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, temp); + xtemp.Add(temp, temp2); + temp2.CopyTo(xtemp); + + gtemp.Multiply(alpha, gtemp); + xtemp.Add(gtemp, temp2); + temp2.CopyTo(xtemp); // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) @@ -498,13 +509,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.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) - zd.Add(d[s].Multiply(beta), zd); + d[s].Multiply(beta, temp); + zd.Add(temp, temp2); + temp2.CopyTo(zd); // z_g = z_g + beta^(jk+i)_((j-1)k+s) g_((j-1)k+s) - zg.Add(g[s].Multiply(beta), zg); + g[s].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); // z_w = z_w + beta^(jk+i)_((j-1)k+s) w_((j-1)k+s) - zw.Add(w[s].Multiply(beta), zw); + w[s].Multiply(beta, temp); + zw.Add(temp, temp2); + temp2.CopyTo(zw); } } @@ -515,14 +532,19 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.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)) - residuals.Add(zw.Multiply(rho), temp); + zw.Multiply(rho, temp2); + residuals.Add(temp2, temp); beta = -_startingVectors[0].DotProduct(temp) / beta; // z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k) - zg.Add(g[k - 1].Multiply(beta), zg); + g[k - 1].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); // z_w = rho_(j+1) (z_w + beta^(jk+i)_((j-1)k+k) w_((j-1)k+k)) - zw.Add(w[k - 1].Multiply(beta), zw); + w[k - 1].Multiply(beta, temp); + zw.Add(temp, temp2); + temp2.CopyTo(zw); zw.Multiply(rho, zw); // z_d = r_(jk+i) + z_w @@ -535,10 +557,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative beta = -_startingVectors[s + 1].DotProduct(zd) / c[s]; // z_d = z_d + beta^(jk+i)_(jk+s) * d_(jk+s) - zd.Add(d[s].Multiply(beta), zd); + d[s].Multiply(beta, temp); + zd.Add(temp, temp2); + temp2.CopyTo(zd); // z_g = z_g + beta^(jk+i)_(jk+s) * g_(jk+s) - zg.Add(g[s].Multiply(beta), zg); + g[s].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); } // d_(jk+i) = z_d - u_(jk+i) @@ -561,19 +587,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative alpha = _startingVectors[i + 1].DotProduct(u) / c[i]; // u_(jk+i+1) = u_(jk+i) - alpha_(jk+i+1) d_(jk+i) - u.Add(d[i].Multiply(-alpha), u); + d[i].Multiply(-alpha, temp); + u.Add(temp, temp2); + temp2.CopyTo(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) - xtemp.Add(gtemp.Multiply(rho * alpha), xtemp); + gtemp.Multiply(rho * alpha, temp); + xtemp.Add(temp, temp2); + temp2.CopyTo(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) - residuals.Add(w[i].Multiply(-rho * alpha), residuals); + w[i].Multiply(-rho * alpha, temp); + residuals.Add(temp, temp2); + temp2.CopyTo(residuals); // We can check the residuals here if they're close if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs index 50bf4fd2..1a2ca2cf 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs @@ -248,6 +248,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -277,7 +282,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative // Temp vectors var temp = new DenseVector(input.Count); - + var temp1 = new DenseVector(input.Count); + var temp2 = new DenseVector(input.Count); + // Initialize var startNorm = input.Norm(2); @@ -319,7 +326,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative alpha = rho / sigma; // yOdd = yEven - alpha * v - yeven.Add(v.Multiply(-alpha), yodd); + v.Multiply(-alpha, temp1); + yeven.Add(temp1, yodd); // Solve M temp = yOdd _preconditioner.Approximate(yodd, temp); @@ -335,7 +343,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative var yinternal = IsEven(iterationNumber) ? yeven : yodd; // pseudoResiduals = pseudoResiduals - alpha * uOdd - pseudoResiduals.Add(uinternal.Multiply(-alpha), pseudoResiduals); + uinternal.Multiply(-alpha, temp1); + pseudoResiduals.Add(temp1, temp2); + temp2.CopyTo(pseudoResiduals); // d = yOdd + theta * theta * eta / alpha * d d.Multiply(theta * theta * eta / alpha, temp); @@ -352,7 +362,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative eta = c * c * alpha; // x = x + eta * d - x.Add(d.Multiply(eta), x); + d.Multiply(eta, temp1); + x.Add(temp1, temp2); + temp2.CopyTo(x); // Check convergence and see if we can bail if (!ShouldContinue(iterationNumber, result, input, pseudoResiduals)) @@ -390,7 +402,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative rho = rhoNew; // yOdd = pseudoResiduals + beta * yOdd - pseudoResiduals.Add(yodd.Multiply(beta), yeven); + yodd.Multiply(beta, temp1); + pseudoResiduals.Add(temp1, yeven); // Solve M temp = yOdd _preconditioner.Approximate(yeven, temp); @@ -399,8 +412,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative matrix.Multiply(temp, ueven); // v = uEven + beta * (uOdd + beta * v) - uodd.Add(v.Multiply(beta), temp); - ueven.Add(temp.Multiply(beta), v); + v.Multiply(beta, temp1); + uodd.Add(temp1, temp); + + temp.Multiply(beta, temp1); + ueven.Add(temp1, v); } // Calculate the real values diff --git a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs index 3eaf2538..ad5862ca 100644 --- a/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Complex32/SparseMatrix.cs @@ -1282,6 +1282,31 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 } } + /// + /// Multiplies this matrix with a vector and places the results into the result vector. + /// + /// The vector to multiply with. + /// The result of the multiplication. + protected override void DoMultiply(Vector rightSide, Vector result) + { + for (var row = 0; row < RowCount; row++) + { + // Get the begin / end index for the current row + var startIndex = _rowIndex[row]; + var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + if (startIndex == endIndex) + { + continue; + } + + var sum = CommonParallel.Aggregate( + startIndex, + endIndex, + index => _nonZeroValues[index] * rightSide[_columnIndices[index]]); + result[row] = sum; + } + } + /// /// Multiplies this matrix with transpose of another matrix and places the results into the result matrix. /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Vector.cs b/src/Numerics/LinearAlgebra/Complex32/Vector.cs index 8c0bf51c..ea0f3dc7 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Vector.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Vector.cs @@ -111,7 +111,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32 /// protected override void DoSubtract(Vector other, Vector result) { - CopyTo(result); CommonParallel.For( 0, Count, diff --git a/src/Numerics/LinearAlgebra/Double/Matrix.cs b/src/Numerics/LinearAlgebra/Double/Matrix.cs index e28d42bf..a79b2998 100644 --- a/src/Numerics/LinearAlgebra/Double/Matrix.cs +++ b/src/Numerics/LinearAlgebra/Double/Matrix.cs @@ -187,7 +187,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double }); } - /// + /// /// Multiplies this matrix with a vector and places the results into the result vector. /// /// The vector to multiply with. diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs index 29310f65..cfc20bb5 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/BiCgStab.cs @@ -258,7 +258,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - + + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -281,9 +286,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative // Choose r~ (for example, r~ = r_0) var tempResiduals = residuals.Clone(); - var temp = result.Clone(); - - // create five temporary vectors needed to hold temporary + + // create seven temporary vectors needed to hold temporary // coefficients. All vectors are mangled in each iteration. // These are defined here to prevent stressing the garbage collector Vector vecP = new DenseVector(residuals.Count); @@ -291,7 +295,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative Vector nu = new DenseVector(residuals.Count); Vector vecS = new DenseVector(residuals.Count); Vector vecSdash = new DenseVector(residuals.Count); - Vector t = new DenseVector(residuals.Count); + Vector temp = new DenseVector(residuals.Count); + Vector temp2 = new DenseVector(residuals.Count); // create some temporary double variables that are needed // to hold values in between iterations @@ -320,10 +325,13 @@ 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)) - vecP.Add(nu.Multiply(-omega), vecP); + nu.Multiply(-omega, temp); + vecP.Add(temp, temp2); + temp2.CopyTo(vecP); vecP.Multiply(beta, vecP); - vecP.Add(residuals, vecP); + vecP.Add(residuals, temp2); + temp2.CopyTo(vecP); } else { @@ -341,7 +349,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative alpha = currentRho * 1 / tempResiduals.DotProduct(nu); // s = r_(i-1) - alpha_i nu_i - residuals.Add(nu.Multiply(-alpha), vecS); + nu.Multiply(-alpha, temp); + residuals.Add(temp, vecS); // Check if we're converged. If so then stop. Otherwise continue; // Calculate the temporary result. @@ -349,8 +358,10 @@ 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); - temp.Add(result, temp); + temp.Add(vecSdash, temp2); + temp2.CopyTo(temp); + temp.Add(result, temp2); + temp2.CopyTo(temp); // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, temp, input, vecS)) @@ -376,18 +387,24 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative _preconditioner.Approximate(vecS, vecSdash); // temp = As~ - matrix.Multiply(vecSdash, t); + matrix.Multiply(vecSdash, temp); // omega_i = temp^T s / temp^T temp - omega = t.DotProduct(vecS) / t.DotProduct(t); + omega = temp.DotProduct(vecS) / temp.DotProduct(temp); // x_i = x_(i-1) + alpha_i p^ + omega_i s^ - result.Add(vecSdash.Multiply(omega), result); - result.Add(vecPdash.Multiply(alpha), result); + temp.Multiply(-omega, residuals); + residuals.Add(vecS, temp2); + temp2.CopyTo(residuals); - t.Multiply(-omega, residuals); - residuals.Add(vecS, residuals); + vecSdash.Multiply(omega, temp); + result.Add(temp, temp2); + temp2.CopyTo(result); + vecPdash.Multiply(alpha, temp); + result.Add(temp, temp2); + temp2.CopyTo(result); + // for continuation it is necessary that omega_i != 0.0 // If omega is only 1 ULP from zero then we fail. if (omega.AlmostEqual(0, 1)) diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs index 07c1a0c4..13bb7ba1 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/GpBiCg.cs @@ -314,8 +314,12 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - // Initialize the solver fields + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) { @@ -362,15 +366,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative Vector z = new DenseVector(residuals.Count); Vector temp = new DenseVector(residuals.Count); - + Vector temp2 = new DenseVector(residuals.Count); + Vector temp3 = new DenseVector(residuals.Count); + // for (k = 0, 1, .... ) var iterationNumber = 0; while (ShouldContinue(iterationNumber, xtemp, input, residuals)) { // p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1)) p.Subtract(u, temp); - - residuals.Add(temp.Multiply(beta), p); + + temp.Multiply(beta, temp2); + residuals.Add(temp2, p); // Solve M b_k = p_k _preconditioner.Approximate(p, temp); @@ -384,15 +391,18 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative // y_k = t_(k-1) - r_k - alpha_k * w_(k-1) + alpha_k s_k s.Subtract(w, temp); t.Subtract(residuals, y); - - y.Add(temp.Multiply(alpha), y); + + temp.Multiply(alpha, temp2); + y.Add(temp2, temp3); + temp3.CopyTo(y); // Store the old value of t in t0 t.CopyTo(t0); // t_k = r_k - alpha_k s_k - residuals.Add(s.Multiply(-alpha), t); - + s.Multiply(-alpha, temp2); + residuals.Add(temp2, t); + // Solve M d_k = t_k _preconditioner.Approximate(t, temp); @@ -449,39 +459,53 @@ 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)) - t0.Add(u.Multiply(beta), temp); + u.Multiply(beta, temp2); + t0.Add(temp2, temp); - temp.Subtract(residuals, temp); + temp.Subtract(residuals, temp3); + temp3.CopyTo(temp); temp.Multiply(eta, temp); - temp.Add(s.Multiply(sigma), u); + s.Multiply(sigma, temp2); + temp.Add(temp2, u); // 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); + u.Multiply(-alpha, temp2); + z.Add(temp2, temp3); + temp3.CopyTo(z); - z.Add(residuals.Multiply(sigma), z); + residuals.Multiply(sigma, temp2); + z.Add(temp2, temp3); + temp3.CopyTo(z); // x_(k+1) = x_k + alpha_k p_k + z_k - xtemp.Add(p.Multiply(alpha), xtemp); - - xtemp.Add(z, xtemp); + p.Multiply(alpha, temp2); + xtemp.Add(temp2, temp3); + temp3.CopyTo(xtemp); + + xtemp.Add(z, temp3); + temp3.CopyTo(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); - t.Add(y.Multiply(-eta), residuals); + y.Multiply(-eta, temp2); + t.Add(temp2, residuals); - residuals.Add(c.Multiply(-sigma), residuals); + c.Multiply(-sigma, temp2); + residuals.Add(temp2, temp3); + temp3.CopyTo(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 - c.Add(s.Multiply(beta), w); + s.Multiply(beta, temp2); + c.Add(temp2, w); // Get the real value _preconditioner.Approximate(xtemp, result); diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs index bf600e04..d4920723 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/MlkBiCgStab.cs @@ -335,6 +335,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -393,6 +398,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative Vector utemp = new DenseVector(residuals.Count); Vector temp = new DenseVector(residuals.Count); Vector temp1 = new DenseVector(residuals.Count); + Vector temp2 = new DenseVector(residuals.Count); Vector zd = new DenseVector(residuals.Count); Vector zg = new DenseVector(residuals.Count); @@ -427,7 +433,8 @@ 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) - residuals.Add(w[k - 1].Multiply(-alpha), u); + w[k - 1].Multiply(-alpha, temp); + residuals.Add(temp, u); // SOLVE M u~_(jk+1) = u_(jk+1) _preconditioner.Approximate(u, temp1); @@ -447,19 +454,23 @@ 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) - xtemp.Add(utemp.Multiply(-rho), xtemp); - - 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, temp); - residuals.Add(temp, residuals); + residuals.Add(temp, temp2); + temp2.CopyTo(residuals); + + // 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, temp); + xtemp.Add(temp, temp2); + temp2.CopyTo(xtemp); + gtemp.Multiply(alpha, gtemp); + xtemp.Add(gtemp, temp2); + temp2.CopyTo(xtemp); + // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) { @@ -497,13 +508,19 @@ 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) - zd.Add(d[s].Multiply(beta), zd); + d[s].Multiply(beta, temp); + zd.Add(temp, temp2); + temp2.CopyTo(zd); // z_g = z_g + beta^(jk+i)_((j-1)k+s) g_((j-1)k+s) - zg.Add(g[s].Multiply(beta), zg); + g[s].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); // z_w = z_w + beta^(jk+i)_((j-1)k+s) w_((j-1)k+s) - zw.Add(w[s].Multiply(beta), zw); + w[s].Multiply(beta, temp); + zw.Add(temp, temp2); + temp2.CopyTo(zw); } } @@ -514,14 +531,19 @@ 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)) - residuals.Add(zw.Multiply(rho), temp); + zw.Multiply(rho, temp2); + residuals.Add(temp2, temp); beta = -_startingVectors[0].DotProduct(temp) / beta; // z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k) - zg.Add(g[k - 1].Multiply(beta), zg); + g[k - 1].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); // z_w = rho_(j+1) (z_w + beta^(jk+i)_((j-1)k+k) w_((j-1)k+k)) - zw.Add(w[k - 1].Multiply(beta), zw); + w[k - 1].Multiply(beta, temp); + zw.Add(temp, temp2); + temp2.CopyTo(zw); zw.Multiply(rho, zw); // z_d = r_(jk+i) + z_w @@ -534,10 +556,14 @@ 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) - zd.Add(d[s].Multiply(beta), zd); + d[s].Multiply(beta, temp); + zd.Add(temp, temp2); + temp2.CopyTo(zd); // z_g = z_g + beta^(jk+i)_(jk+s) * g_(jk+s) - zg.Add(g[s].Multiply(beta), zg); + g[s].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); } // d_(jk+i) = z_d - u_(jk+i) @@ -560,19 +586,25 @@ 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) - u.Add(d[i].Multiply(-alpha), u); + d[i].Multiply(-alpha, temp); + u.Add(temp, temp2); + temp2.CopyTo(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) - xtemp.Add(gtemp.Multiply(rho * alpha), xtemp); + gtemp.Multiply(rho * alpha, temp); + xtemp.Add(temp, temp2); + temp2.CopyTo(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) - residuals.Add(w[i].Multiply(-rho * alpha), residuals); + w[i].Multiply(-rho * alpha, temp); + residuals.Add(temp, temp2); + temp2.CopyTo(residuals); // We can check the residuals here if they're close if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs index 5720ce95..2b9d1c38 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs @@ -247,6 +247,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -276,6 +281,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative // Temp vectors var temp = new DenseVector(input.Count); + var temp1 = new DenseVector(input.Count); + var temp2 = new DenseVector(input.Count); // Initialize var startNorm = input.Norm(2); @@ -318,7 +325,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative alpha = rho / sigma; // yOdd = yEven - alpha * v - yeven.Add(v.Multiply(-alpha), yodd); + v.Multiply(-alpha, temp1); + yeven.Add(temp1, yodd); // Solve M temp = yOdd _preconditioner.Approximate(yodd, temp); @@ -334,7 +342,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative var yinternal = IsEven(iterationNumber) ? yeven : yodd; // pseudoResiduals = pseudoResiduals - alpha * uOdd - pseudoResiduals.Add(uinternal.Multiply(-alpha), pseudoResiduals); + uinternal.Multiply(-alpha, temp1); + pseudoResiduals.Add(temp1, temp2); + temp2.CopyTo(pseudoResiduals); // d = yOdd + theta * theta * eta / alpha * d d.Multiply(theta * theta * eta / alpha, temp); @@ -351,7 +361,9 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative eta = c * c * alpha; // x = x + eta * d - x.Add(d.Multiply(eta), x); + d.Multiply(eta, temp1); + x.Add(temp1, temp2); + temp2.CopyTo(x); // Check convergence and see if we can bail if (!ShouldContinue(iterationNumber, result, input, pseudoResiduals)) @@ -389,7 +401,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative rho = rhoNew; // yOdd = pseudoResiduals + beta * yOdd - pseudoResiduals.Add(yodd.Multiply(beta), yeven); + yodd.Multiply(beta, temp1); + pseudoResiduals.Add(temp1, yeven); // Solve M temp = yOdd _preconditioner.Approximate(yeven, temp); @@ -398,9 +411,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative matrix.Multiply(temp, ueven); // v = uEven + beta * (uOdd + beta * v) - uodd.Add(v.Multiply(beta), temp); + v.Multiply(beta, temp1); + uodd.Add(temp1, temp); - ueven.Add(temp.Multiply(beta), v); + temp.Multiply(beta, temp1); + ueven.Add(temp1, v); } // Calculate the real values diff --git a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs index 61e832de..8207a234 100644 --- a/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Double/SparseMatrix.cs @@ -1281,6 +1281,31 @@ namespace MathNet.Numerics.LinearAlgebra.Double } } + /// + /// Multiplies this matrix with a vector and places the results into the result vector. + /// + /// The vector to multiply with. + /// The result of the multiplication. + protected override void DoMultiply(Vector rightSide, Vector result) + { + for (var row = 0; row < RowCount; row++) + { + // Get the begin / end index for the current row + var startIndex = _rowIndex[row]; + var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + if (startIndex == endIndex) + { + continue; + } + + var sum = CommonParallel.Aggregate( + startIndex, + endIndex, + index => _nonZeroValues[index] * rightSide[_columnIndices[index]]); + result[row] = sum; + } + } + /// /// Multiplies this matrix with transpose of another matrix and places the results into the result matrix. /// diff --git a/src/Numerics/LinearAlgebra/Double/Vector.cs b/src/Numerics/LinearAlgebra/Double/Vector.cs index 1de2b78a..7baefb6a 100644 --- a/src/Numerics/LinearAlgebra/Double/Vector.cs +++ b/src/Numerics/LinearAlgebra/Double/Vector.cs @@ -110,7 +110,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double /// protected override void DoSubtract(Vector other, Vector result) { - CopyTo(result); CommonParallel.For( 0, Count, diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/BiCgStab.cs b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/BiCgStab.cs index b37a286a..3b2c7a63 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/BiCgStab.cs @@ -258,7 +258,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative { throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - + + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -281,9 +286,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative // Choose r~ (for example, r~ = r_0) var tempResiduals = residuals.Clone(); - var temp = result.Clone(); - // create five temporary vectors needed to hold temporary + // create seven temporary vectors needed to hold temporary // coefficients. All vectors are mangled in each iteration. // These are defined here to prevent stressing the garbage collector Vector vecP = new DenseVector(residuals.Count); @@ -291,7 +295,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative Vector nu = new DenseVector(residuals.Count); Vector vecS = new DenseVector(residuals.Count); Vector vecSdash = new DenseVector(residuals.Count); - Vector t = new DenseVector(residuals.Count); + Vector temp = new DenseVector(residuals.Count); + Vector temp2 = new DenseVector(residuals.Count); // create some temporary float variables that are needed // to hold values in between iterations @@ -320,10 +325,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.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)) - vecP.Add(nu.Multiply(-omega), vecP); + nu.Multiply(-omega, temp); + vecP.Add(temp, temp2); + temp2.CopyTo(vecP); vecP.Multiply(beta, vecP); - vecP.Add(residuals, vecP); + vecP.Add(residuals, temp2); + temp2.CopyTo(vecP); } else { @@ -341,7 +349,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative alpha = currentRho * 1 / tempResiduals.DotProduct(nu); // s = r_(i-1) - alpha_i nu_i - residuals.Add(nu.Multiply(-alpha), vecS); + nu.Multiply(-alpha, temp); + residuals.Add(temp, vecS); // Check if we're converged. If so then stop. Otherwise continue; // Calculate the temporary result. @@ -349,8 +358,10 @@ namespace MathNet.Numerics.LinearAlgebra.Single.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); - temp.Add(result, temp); + temp.Add(vecSdash, temp2); + temp2.CopyTo(temp); + temp.Add(result, temp2); + temp2.CopyTo(temp); // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, temp, input, vecS)) @@ -376,17 +387,23 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative _preconditioner.Approximate(vecS, vecSdash); // temp = As~ - matrix.Multiply(vecSdash, t); + matrix.Multiply(vecSdash, temp); // omega_i = temp^T s / temp^T temp - omega = t.DotProduct(vecS) / t.DotProduct(t); + omega = temp.DotProduct(vecS) / temp.DotProduct(temp); // x_i = x_(i-1) + alpha_i p^ + omega_i s^ - result.Add(vecSdash.Multiply(omega), result); - result.Add(vecPdash.Multiply(alpha), result); + temp.Multiply(-omega, residuals); + residuals.Add(vecS, temp2); + temp2.CopyTo(residuals); + + vecSdash.Multiply(omega, temp); + result.Add(temp, temp2); + temp2.CopyTo(result); - t.Multiply(-omega, residuals); - residuals.Add(vecS, residuals); + vecPdash.Multiply(alpha, temp); + result.Add(temp, temp2); + temp2.CopyTo(result); // for continuation it is necessary that omega_i != 0.0 // If omega is only 1 ULP from zero then we fail. diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/GpBiCg.cs b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/GpBiCg.cs index 10de0666..0e8e52ed 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/GpBiCg.cs @@ -314,8 +314,12 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } - // Initialize the solver fields + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) { @@ -362,15 +366,18 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative Vector z = new DenseVector(residuals.Count); Vector temp = new DenseVector(residuals.Count); - + Vector temp2 = new DenseVector(residuals.Count); + Vector temp3 = new DenseVector(residuals.Count); + // for (k = 0, 1, .... ) var iterationNumber = 0; while (ShouldContinue(iterationNumber, xtemp, input, residuals)) { // p_k = r_k + beta_(k-1) * (p_(k-1) - u_(k-1)) p.Subtract(u, temp); - - residuals.Add(temp.Multiply(beta), p); + + temp.Multiply(beta, temp2); + residuals.Add(temp2, p); // Solve M b_k = p_k _preconditioner.Approximate(p, temp); @@ -384,14 +391,17 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative // y_k = t_(k-1) - r_k - alpha_k * w_(k-1) + alpha_k s_k s.Subtract(w, temp); t.Subtract(residuals, y); - - y.Add(temp.Multiply(alpha), y); + + temp.Multiply(alpha, temp2); + y.Add(temp2, temp3); + temp3.CopyTo(y); // Store the old value of t in t0 t.CopyTo(t0); // t_k = r_k - alpha_k s_k - residuals.Add(s.Multiply(-alpha), t); + s.Multiply(-alpha, temp2); + residuals.Add(temp2, t); // Solve M d_k = t_k _preconditioner.Approximate(t, temp); @@ -449,39 +459,53 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative } // u_k = sigma_k s_k + eta_k (t_(k-1) - r_k + beta_(k-1) u_(k-1)) - t0.Add(u.Multiply(beta), temp); + u.Multiply(beta, temp2); + t0.Add(temp2, temp); - temp.Subtract(residuals, temp); + temp.Subtract(residuals, temp3); + temp3.CopyTo(temp); temp.Multiply(eta, temp); - temp.Add(s.Multiply(sigma), u); + s.Multiply(sigma, temp2); + temp.Add(temp2, u); // 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); + u.Multiply(-alpha, temp2); + z.Add(temp2, temp3); + temp3.CopyTo(z); - z.Add(residuals.Multiply(sigma), z); + residuals.Multiply(sigma, temp2); + z.Add(temp2, temp3); + temp3.CopyTo(z); // x_(k+1) = x_k + alpha_k p_k + z_k - xtemp.Add(p.Multiply(alpha), xtemp); + p.Multiply(alpha, temp2); + xtemp.Add(temp2, temp3); + temp3.CopyTo(xtemp); - xtemp.Add(z, xtemp); + xtemp.Add(z, temp3); + temp3.CopyTo(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); - t.Add(y.Multiply(-eta), residuals); + y.Multiply(-eta, temp2); + t.Add(temp2, residuals); - residuals.Add(c.Multiply(-sigma), residuals); + c.Multiply(-sigma, temp2); + residuals.Add(temp2, temp3); + temp3.CopyTo(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 - c.Add(s.Multiply(beta), w); + s.Multiply(beta, temp2); + c.Add(temp2, w); // Get the real value _preconditioner.Approximate(xtemp, result); diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/MlkBiCgStab.cs index 2b959c28..08688999 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/MlkBiCgStab.cs @@ -334,6 +334,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -392,6 +397,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative Vector utemp = new DenseVector(residuals.Count); Vector temp = new DenseVector(residuals.Count); Vector temp1 = new DenseVector(residuals.Count); + Vector temp2 = new DenseVector(residuals.Count); Vector zd = new DenseVector(residuals.Count); Vector zg = new DenseVector(residuals.Count); @@ -426,7 +432,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.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) - residuals.Add(w[k - 1].Multiply(-alpha), u); + w[k - 1].Multiply(-alpha, temp); + residuals.Add(temp, u); // SOLVE M u~_(jk+1) = u_(jk+1) _preconditioner.Approximate(u, temp1); @@ -446,18 +453,22 @@ namespace MathNet.Numerics.LinearAlgebra.Single.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) - xtemp.Add(utemp.Multiply(-rho), xtemp); - - 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, temp); - residuals.Add(temp, residuals); + residuals.Add(temp, temp2); + temp2.CopyTo(residuals); + + // 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, temp); + xtemp.Add(temp, temp2); + temp2.CopyTo(xtemp); + + gtemp.Multiply(alpha, gtemp); + xtemp.Add(gtemp, temp2); + temp2.CopyTo(xtemp); // Check convergence and stop if we are converged. if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) @@ -496,13 +507,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.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) - zd.Add(d[s].Multiply(beta), zd); + d[s].Multiply(beta, temp); + zd.Add(temp, temp2); + temp2.CopyTo(zd); // z_g = z_g + beta^(jk+i)_((j-1)k+s) g_((j-1)k+s) - zg.Add(g[s].Multiply(beta), zg); + g[s].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); // z_w = z_w + beta^(jk+i)_((j-1)k+s) w_((j-1)k+s) - zw.Add(w[s].Multiply(beta), zw); + w[s].Multiply(beta, temp); + zw.Add(temp, temp2); + temp2.CopyTo(zw); } } @@ -513,14 +530,19 @@ namespace MathNet.Numerics.LinearAlgebra.Single.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)) - residuals.Add(zw.Multiply(rho), temp); + zw.Multiply(rho, temp2); + residuals.Add(temp2, temp); beta = -_startingVectors[0].DotProduct(temp) / beta; // z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k) - zg.Add(g[k - 1].Multiply(beta), zg); + g[k - 1].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); // z_w = rho_(j+1) (z_w + beta^(jk+i)_((j-1)k+k) w_((j-1)k+k)) - zw.Add(w[k - 1].Multiply(beta), zw); + w[k - 1].Multiply(beta, temp); + zw.Add(temp, temp2); + temp2.CopyTo(zw); zw.Multiply(rho, zw); // z_d = r_(jk+i) + z_w @@ -533,10 +555,14 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative beta = -_startingVectors[s + 1].DotProduct(zd) / c[s]; // z_d = z_d + beta^(jk+i)_(jk+s) * d_(jk+s) - zd.Add(d[s].Multiply(beta), zd); + d[s].Multiply(beta, temp); + zd.Add(temp, temp2); + temp2.CopyTo(zd); // z_g = z_g + beta^(jk+i)_(jk+s) * g_(jk+s) - zg.Add(g[s].Multiply(beta), zg); + g[s].Multiply(beta, temp); + zg.Add(temp, temp2); + temp2.CopyTo(zg); } // d_(jk+i) = z_d - u_(jk+i) @@ -559,19 +585,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative alpha = _startingVectors[i + 1].DotProduct(u) / c[i]; // u_(jk+i+1) = u_(jk+i) - alpha_(jk+i+1) d_(jk+i) - u.Add(d[i].Multiply(-alpha), u); + d[i].Multiply(-alpha, temp); + u.Add(temp, temp2); + temp2.CopyTo(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) - xtemp.Add(gtemp.Multiply(rho * alpha), xtemp); + gtemp.Multiply(rho * alpha, temp); + xtemp.Add(temp, temp2); + temp2.CopyTo(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) - residuals.Add(w[i].Multiply(-rho * alpha), residuals); + w[i].Multiply(-rho * alpha, temp); + residuals.Add(temp, temp2); + temp2.CopyTo(residuals); // We can check the residuals here if they're close if (!ShouldContinue(iterationNumber, xtemp, input, residuals)) diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs index 0f388ea0..f5191cc5 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs @@ -247,6 +247,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative throw new ArgumentException(Resources.ArgumentVectorsSameLength); } + if (input.Count != matrix.RowCount) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + // Initialize the solver fields // Set the convergence monitor if (_iterator == null) @@ -276,7 +281,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative // Temp vectors var temp = new DenseVector(input.Count); - + var temp1 = new DenseVector(input.Count); + var temp2 = new DenseVector(input.Count); + // Initialize var startNorm = (float)input.Norm(2); @@ -318,7 +325,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative alpha = rho / sigma; // yOdd = yEven - alpha * v - yeven.Add(v.Multiply(-alpha), yodd); + v.Multiply(-alpha, temp1); + yeven.Add(temp1, yodd); // Solve M temp = yOdd _preconditioner.Approximate(yodd, temp); @@ -334,7 +342,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative var yinternal = IsEven(iterationNumber) ? yeven : yodd; // pseudoResiduals = pseudoResiduals - alpha * uOdd - pseudoResiduals.Add(uinternal.Multiply(-alpha), pseudoResiduals); + uinternal.Multiply(-alpha, temp1); + pseudoResiduals.Add(temp1, temp2); + temp2.CopyTo(pseudoResiduals); // d = yOdd + theta * theta * eta / alpha * d d.Multiply(theta * theta * eta / alpha, temp); @@ -351,7 +361,9 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative eta = c * c * alpha; // x = x + eta * d - x.Add(d.Multiply(eta), x); + d.Multiply(eta, temp1); + x.Add(temp1, temp2); + temp2.CopyTo(x); // Check convergence and see if we can bail if (!ShouldContinue(iterationNumber, result, input, pseudoResiduals)) @@ -389,7 +401,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative rho = rhoNew; // yOdd = pseudoResiduals + beta * yOdd - pseudoResiduals.Add(yodd.Multiply(beta), yeven); + yodd.Multiply(beta, temp1); + pseudoResiduals.Add(temp1, yeven); // Solve M temp = yOdd _preconditioner.Approximate(yeven, temp); @@ -398,9 +411,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative matrix.Multiply(temp, ueven); // v = uEven + beta * (uOdd + beta * v) - uodd.Add(v.Multiply(beta), temp); + v.Multiply(beta, temp1); + uodd.Add(temp1, temp); - ueven.Add(temp.Multiply(beta), v); + temp.Multiply(beta, temp1); + ueven.Add(temp1, v); } // Calculate the real values diff --git a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs index 0f81c913..d79efe6c 100644 --- a/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs +++ b/src/Numerics/LinearAlgebra/Single/SparseMatrix.cs @@ -1281,6 +1281,31 @@ namespace MathNet.Numerics.LinearAlgebra.Single } } + /// + /// Multiplies this matrix with a vector and places the results into the result vector. + /// + /// The vector to multiply with. + /// The result of the multiplication. + protected override void DoMultiply(Vector rightSide, Vector result) + { + for (var row = 0; row < RowCount; row++) + { + // Get the begin / end index for the current row + var startIndex = _rowIndex[row]; + var endIndex = row < _rowIndex.Length - 1 ? _rowIndex[row + 1] : NonZerosCount; + if (startIndex == endIndex) + { + continue; + } + + var sum = CommonParallel.Aggregate( + startIndex, + endIndex, + index => _nonZeroValues[index] * rightSide[_columnIndices[index]]); + result[row] = sum; + } + } + /// /// Multiplies this matrix with transpose of another matrix and places the results into the result matrix. /// diff --git a/src/Numerics/LinearAlgebra/Single/Vector.cs b/src/Numerics/LinearAlgebra/Single/Vector.cs index fe6c1378..0248c1ac 100644 --- a/src/Numerics/LinearAlgebra/Single/Vector.cs +++ b/src/Numerics/LinearAlgebra/Single/Vector.cs @@ -111,7 +111,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single /// protected override void DoSubtract(Vector other, Vector result) { - CopyTo(result); CommonParallel.For( 0, Count, @@ -467,7 +466,5 @@ namespace MathNet.Numerics.LinearAlgebra.Single return v; } - - } } diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/BiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/BiCgStabTest.cs index 3b50131b..c1a22c6b 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/BiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/BiCgStabTest.cs @@ -86,7 +86,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Solvers.Iterativ Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply(new Complex(System.Math.PI, System.Math.PI)); + matrix.Multiply(new Complex(System.Math.PI, System.Math.PI), matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, Complex.One); diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/GpBiCgTest.cs b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/GpBiCgTest.cs index 97de96f7..63926d5a 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/GpBiCgTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/GpBiCgTest.cs @@ -86,7 +86,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Solvers.Iterativ Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply(System.Math.PI); + matrix.Multiply(System.Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/MlkBiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/MlkBiCgStabTest.cs index 290fdc92..0a10c267 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/MlkBiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/MlkBiCgStabTest.cs @@ -87,7 +87,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Solvers.Iterativ Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply(System.Math.PI); + matrix.Multiply(System.Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/TFQMRTest.cs b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/TFQMRTest.cs index 532b63bd..7bc5f74a 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/TFQMRTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex/Solvers/Iterative/TFQMRTest.cs @@ -87,7 +87,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex.Solvers.Iterativ Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply(System.Math.PI); + matrix.Multiply(System.Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/BiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/BiCgStabTest.cs index 3fd13dc4..5017436c 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/BiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/BiCgStabTest.cs @@ -87,7 +87,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.Iterat Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply(new Complex32((float)Math.PI, (float)Math.PI)); + matrix.Multiply(new Complex32((float)Math.PI, (float)Math.PI), matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, Complex32.One); diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/GpBiCgTest.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/GpBiCgTest.cs index 9011d14f..93185795 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/GpBiCgTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/GpBiCgTest.cs @@ -87,7 +87,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.Iterat Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply((float)Math.PI); + matrix.Multiply((float)Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/MlkBiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/MlkBiCgStabTest.cs index c1aca5ba..a73534b5 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/MlkBiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/MlkBiCgStabTest.cs @@ -88,7 +88,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.Iterat Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply((float)Math.PI); + matrix.Multiply((float)Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/TFQMRTest.cs b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/TFQMRTest.cs index d1858690..39b65ce9 100644 --- a/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/TFQMRTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Complex32/Solvers/Iterative/TFQMRTest.cs @@ -88,7 +88,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Complex32.Solvers.Iterat Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply((float)Math.PI); + matrix.Multiply((float)Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs index daa9a26f..96bcd428 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/BiCgStabTest.cs @@ -85,7 +85,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply(System.Math.PI); + matrix.Multiply(System.Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs index 666584d7..96f706ce 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/GpBiCgTest.cs @@ -85,7 +85,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply(System.Math.PI); + matrix.Multiply(System.Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs index 188fb871..efc9a5f9 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/MlkBiCgStabTest.cs @@ -86,7 +86,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply(System.Math.PI); + matrix.Multiply(System.Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs index 17b5b5cf..3919f37d 100644 --- a/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Double/Solvers/Iterative/TFQMRTest.cs @@ -86,7 +86,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Solvers.Iterative Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply(System.Math.PI); + matrix.Multiply(System.Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/BiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/BiCgStabTest.cs index 9649eb4e..1f97b779 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/BiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/BiCgStabTest.cs @@ -86,7 +86,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.Iterative Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply((float)Math.PI); + matrix.Multiply((float)Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/GpBiCgTest.cs b/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/GpBiCgTest.cs index e49c80d0..77f83524 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/GpBiCgTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/GpBiCgTest.cs @@ -86,7 +86,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.Iterative Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply((float)Math.PI); + matrix.Multiply((float)Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/MlkBiCgStabTest.cs b/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/MlkBiCgStabTest.cs index 829ea5de..98c8c68c 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/MlkBiCgStabTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/MlkBiCgStabTest.cs @@ -87,7 +87,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.Iterative Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply((float)Math.PI); + matrix.Multiply((float)Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1); diff --git a/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/TFQMRTest.cs b/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/TFQMRTest.cs index d0f98012..0825e408 100644 --- a/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/TFQMRTest.cs +++ b/src/UnitTests/LinearAlgebraTests/Single/Solvers/Iterative/TFQMRTest.cs @@ -87,7 +87,7 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Solvers.Iterative Matrix matrix = SparseMatrix.Identity(100); // Scale it with a funny number - matrix.Multiply((float)Math.PI); + matrix.Multiply((float)Math.PI, matrix); // Create the y vector Vector y = new DenseVector(matrix.RowCount, 1);