diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs index 2293b2ee..bdfbb674 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs @@ -39,6 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif /// @@ -76,24 +77,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// The status used if there is no status, i.e. the solver hasn't run yet and there is no /// iterator. /// - private static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); + static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); /// /// The preconditioner that will be used. Can be set to , in which case the default /// pre-conditioner will be used. /// - private IPreConditioner _preconditioner; + IPreConditioner _preconditioner; /// /// The iterative process controller. /// - private IIterator _iterator; + IIterator _iterator; /// /// Indicates if the user has stopped the solver. /// - private bool _hasBeenStopped; - + bool _hasBeenStopped; + /// /// Initializes a new instance of the class. /// @@ -101,7 +102,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// When using this constructor the solver will use the with /// the standard settings and a default preconditioner. /// - public BiCgStab() : this(null, null) + public BiCgStab() + : this(null, null) { } @@ -124,7 +126,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// /// /// The that will be used to monitor the iterative process. - public BiCgStab(IIterator iterator) : this(null, iterator) + public BiCgStab(IIterator iterator) + : this(null, iterator) { } @@ -136,7 +139,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// the standard settings. /// /// The that will be used to precondition the matrix equation. - public BiCgStab(IPreConditioner preconditioner) : this(preconditioner, null) + public BiCgStab(IPreConditioner preconditioner) + : this(preconditioner, null) { } @@ -186,10 +190,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// public ICalculationStatus IterationResult { - get - { - return (_iterator != null) ? _iterator.Status : DefaultStatus; - } + get { return (_iterator != null) ? _iterator.Status : DefaultStatus; } } /// @@ -278,9 +279,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative { _preconditioner = new UnitPreconditioner(); } - + _preconditioner.Initialize(matrix); - + // Compute r_0 = b - Ax_0 for some initial guess x_0 // In this case we take x_0 = vector // This is basically a SAXPY so it could be made a lot faster @@ -312,7 +313,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative { // rho_(i-1) = r~^T r_(i-1) // dotproduct r~ and r_(i-1) var oldRho = currentRho; - currentRho = tempResiduals.DotProduct(residuals); + currentRho = tempResiduals.ConjugateDotProduct(residuals); // if (rho_(i-1) == 0) // METHOD FAILS // If rho is only 1 ULP from zero then we fail. @@ -325,7 +326,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative if (iterationNumber != 0) { // beta_(i-1) = (rho_(i-1)/rho_(i-2))(alpha_(i-1)/omega(i-1)) - var beta = (currentRho / oldRho) * (alpha / omega); + 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, temp); @@ -344,12 +345,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative // SOLVE Mp~ = p_i // M = preconditioner _preconditioner.Approximate(vecP, vecPdash); - + // nu_i = Ap~ matrix.Multiply(vecPdash, nu); // alpha_i = rho_(i-1)/ (r~^T nu_i) = rho / dotproduct(r~ and nu_i) - alpha = currentRho * 1 / tempResiduals.DotProduct(nu); + alpha = currentRho*1/tempResiduals.ConjugateDotProduct(nu); // s = r_(i-1) - alpha_i nu_i nu.Multiply(-alpha, temp); @@ -393,7 +394,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative matrix.Multiply(vecSdash, temp); // omega_i = temp^T s / temp^T temp - omega = temp.DotProduct(vecS) / temp.DotProduct(temp); + omega = temp.ConjugateDotProduct(vecS)/temp.ConjugateDotProduct(temp); // x_i = x_(i-1) + alpha_i p^ + omega_i s^ temp.Multiply(-omega, residuals); @@ -437,11 +438,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// Residual values in . /// Instance of the x. /// Instance of the b. - private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { // -Ax = residual matrix.Multiply(x, residual); - + // Do not use residual = residual.Negate() because it creates another object residual.Multiply(-1, residual); @@ -457,7 +458,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// Source . /// Residual . /// true if continue, otherwise false - private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) + bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) { if (_hasBeenStopped) { @@ -493,7 +494,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative throw new ArgumentNullException("input"); } - var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount); + var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount); Solve(matrix, input, result); return result; } @@ -529,7 +530,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, (Vector)input.Column(column)); + var solution = Solve(matrix, (Vector) input.Column(column)); foreach (var element in solution.GetIndexedEnumerator()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs index 2584d40a..a866c214 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs @@ -39,6 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif /// @@ -377,7 +378,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative matrix.Multiply(temp, s); // alpha_k = (r*_0 * r_k) / (r*_0 * s_k) - var alpha = rdash.DotProduct(residuals)/rdash.DotProduct(s); + var alpha = rdash.ConjugateDotProduct(residuals)/rdash.ConjugateDotProduct(s); // y_k = t_(k-1) - r_k - alpha_k * w_(k-1) + alpha_k s_k s.Subtract(w, temp); @@ -399,7 +400,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative // c_k = A d_k matrix.Multiply(temp, c); - var cdot = c.DotProduct(c); + var cdot = c.ConjugateDotProduct(c); // cDot can only be zero if c is a zero vector // We'll set cDot to 1 if it is zero to prevent NaN's @@ -414,7 +415,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative // to do at least one at the start to initialize the // system, but we'll only have to take special measures // if we don't do any so ... - var ctdot = c.DotProduct(t); + var ctdot = c.ConjugateDotProduct(t); Complex eta; Complex sigma; if (((_numberOfBiCgStabSteps == 0) && (iterationNumber == 0)) || ShouldRunBiCgStabSteps(iterationNumber)) @@ -427,7 +428,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative } else { - var ydot = y.DotProduct(y); + var ydot = y.ConjugateDotProduct(y); // yDot can only be zero if y is a zero vector // We'll set yDot to 1 if it is zero to prevent NaN's @@ -438,8 +439,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative ydot = 1.0; } - var ytdot = y.DotProduct(t); - var cydot = c.DotProduct(y); + var ytdot = y.ConjugateDotProduct(t); + var cydot = c.ConjugateDotProduct(y); var denom = (cdot*ydot) - (cydot*cydot); @@ -493,7 +494,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative // 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; + beta = (!sigma.Real.AlmostEqual(0, 1) || !sigma.Imaginary.AlmostEqual(0, 1)) ? alpha/sigma*rdash.ConjugateDotProduct(residuals)/rdash.ConjugateDotProduct(t0) : 0; // w_k = c_k + beta_k s_k s.Multiply(beta, temp2); diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs index 8e4c9ccc..2f326733 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs @@ -43,6 +43,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif /// @@ -73,39 +74,39 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// /// The default number of starting vectors. /// - private const int DefaultNumberOfStartingVectors = 50; - + const int DefaultNumberOfStartingVectors = 50; + /// /// The status used if there is no status, i.e. the solver hasn't run yet and there is no /// iterator. /// - private static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); + static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); /// /// The preconditioner that will be used. Can be set to , in which case the default /// pre-conditioner will be used. /// - private IPreConditioner _preconditioner; + IPreConditioner _preconditioner; /// /// The iterative process controller. /// - private IIterator _iterator; + IIterator _iterator; /// /// The collection of starting vectors which are used as the basis for the Krylov sub-space. /// - private IList _startingVectors; + IList _startingVectors; /// /// The number of starting vectors used by the algorithm /// - private int _numberOfStartingVectors = DefaultNumberOfStartingVectors; + int _numberOfStartingVectors = DefaultNumberOfStartingVectors; /// /// Indicates if the user has stopped the solver. /// - private bool _hasBeenStopped; + bool _hasBeenStopped; /// /// Initializes a new instance of the class. @@ -114,7 +115,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// When using this constructor the solver will use the with /// the standard settings and a default preconditioner. /// - public MlkBiCgStab() : this(null, null) + public MlkBiCgStab() + : this(null, null) { } @@ -137,7 +139,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// /// /// The that will be used to monitor the iterative process. - public MlkBiCgStab(IIterator iterator) : this(null, iterator) + public MlkBiCgStab(IIterator iterator) + : this(null, iterator) { } @@ -149,7 +152,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// the standard settings. /// /// The that will be used to precondition the matrix equation. - public MlkBiCgStab(IPreConditioner preconditioner) : this(preconditioner, null) + public MlkBiCgStab(IPreConditioner preconditioner) + : this(preconditioner, null) { } @@ -186,10 +190,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative public int NumberOfStartingVectors { [DebuggerStepThrough] - get - { - return _numberOfStartingVectors; - } + get { return _numberOfStartingVectors; } [DebuggerStepThrough] set @@ -236,10 +237,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative public IList StartingVectors { [DebuggerStepThrough] - get - { - return _startingVectors; - } + get { return _startingVectors; } [DebuggerStepThrough] set @@ -261,10 +259,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative public ICalculationStatus IterationResult { [DebuggerStepThrough] - get - { - return (_iterator != null) ? _iterator.Status : DefaultStatus; - } + get { return (_iterator != null) ? _iterator.Status : DefaultStatus; } } /// @@ -348,7 +343,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative { _preconditioner = new UnitPreconditioner(); } - + _preconditioner.Initialize(matrix); // Choose an initial guess x_0 @@ -402,7 +397,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative Vector zw = new DenseVector(residuals.Count); var d = CreateVectorArray(_startingVectors.Count, residuals.Count); - + // g_0 = r_0 var g = CreateVectorArray(_startingVectors.Count, residuals.Count); residuals.CopyTo(g[k - 1]); @@ -420,14 +415,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative matrix.Multiply(gtemp, w[k - 1]); // c_((j-1)k+k) = q^T_1 w_((j-1)k+k) - c[k - 1] = _startingVectors[0].DotProduct(w[k - 1]); + c[k - 1] = _startingVectors[0].ConjugateDotProduct(w[k - 1]); if (c[k - 1].Real.AlmostEqual(0, 1) && c[k - 1].Imaginary.AlmostEqual(0, 1)) { throw new Exception("Iterative solver experience a numerical break down"); } // alpha_(jk+1) = q^T_1 r_((j-1)k+k) / c_((j-1)k+k) - var alpha = _startingVectors[0].DotProduct(residuals) / c[k - 1]; + var alpha = _startingVectors[0].ConjugateDotProduct(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, temp); @@ -439,7 +434,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative // rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2 matrix.Multiply(temp1, temp); - var rho = temp.DotProduct(temp); + var rho = temp.ConjugateDotProduct(temp); // If rho is zero then temp is a zero vector and we're probably // about to have zero residuals (i.e. an exact solution). @@ -449,7 +444,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative rho = 1.0; } - rho = -u.DotProduct(temp) / rho; + rho = -u.ConjugateDotProduct(temp)/rho; // r_(jk+1) = rho_(j+1) A u~_(jk+1) + u_(jk+1) u.CopyTo(residuals); @@ -502,7 +497,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative for (var s = i; s < k - 1; s++) { // beta^(jk+i)_((j-1)k+s) = -q^t_(s+1) z_d / c_((j-1)k+s) - beta = -_startingVectors[s + 1].DotProduct(zd) / c[s]; + beta = -_startingVectors[s + 1].ConjugateDotProduct(zd)/c[s]; // z_d = z_d + beta^(jk+i)_((j-1)k+s) d_((j-1)k+s) d[s].Multiply(beta, temp); @@ -521,7 +516,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative } } - beta = rho * c[k - 1]; + beta = rho*c[k - 1]; if (beta.Real.AlmostEqual(0, 1) && beta.Imaginary.AlmostEqual(0, 1)) { throw new Exception("Iterative solver experience a numerical break down"); @@ -530,7 +525,7 @@ 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)) zw.Multiply(rho, temp2); residuals.Add(temp2, temp); - beta = -_startingVectors[0].DotProduct(temp) / beta; + beta = -_startingVectors[0].ConjugateDotProduct(temp)/beta; // z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k) g[k - 1].Multiply(beta, temp); @@ -550,7 +545,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative for (var s = 0; s < i - 1; s++) { // beta^(jk+i)_(jk+s) = -q^T_s+1 z_d / c_(jk+s) - beta = -_startingVectors[s + 1].DotProduct(zd) / c[s]; + beta = -_startingVectors[s + 1].ConjugateDotProduct(zd)/c[s]; // z_d = z_d + beta^(jk+i)_(jk+s) * d_(jk+s) d[s].Multiply(beta, temp); @@ -573,14 +568,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative if (i < k - 1) { // c_(jk+1) = q^T_i+1 d_(jk+i) - c[i] = _startingVectors[i + 1].DotProduct(d[i]); + c[i] = _startingVectors[i + 1].ConjugateDotProduct(d[i]); if (c[i].Real.AlmostEqual(0, 1) && c[i].Imaginary.AlmostEqual(0, 1)) { throw new Exception("Iterative solver experience a numerical break down"); } // alpha_(jk+i+1) = q^T_(i+1) u_(jk+i) / c_(jk+i) - alpha = _startingVectors[i + 1].DotProduct(u) / c[i]; + alpha = _startingVectors[i + 1].ConjugateDotProduct(u)/c[i]; // u_(jk+i+1) = u_(jk+i) - alpha_(jk+i+1) d_(jk+i) d[i].Multiply(-alpha, temp); @@ -591,7 +586,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative _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, temp); + gtemp.Multiply(rho*alpha, temp); xtemp.Add(temp, temp2); temp2.CopyTo(xtemp); @@ -599,7 +594,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative 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, temp); + w[i].Multiply(-rho*alpha, temp); residuals.Add(temp, temp2); temp2.CopyTo(residuals); @@ -626,7 +621,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// Maximum number /// Number of variables /// Number of starting vectors to create - private static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables) + static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables) { // Create no more starting vectors than the size of the problem - 1 return Math.Min(maximumNumberOfStartingVectors, (numberOfVariables - 1)); @@ -643,7 +638,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// the is smaller than /// the . /// - private static IList CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) + static IList CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) { // Create no more starting vectors than the size of the problem - 1 // Get random values and then orthogonalize them with @@ -677,10 +672,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative var result = new List(); for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) { - result.Add((Vector)orthogonalMatrix.Column(i)); - + result.Add((Vector) orthogonalMatrix.Column(i)); + // Normalize the result vector - result[i].Multiply(1 / result[i].L2Norm(), result[i]); + result[i].Multiply(1/result[i].L2Norm(), result[i]); } return result; @@ -692,7 +687,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// Number of vectors /// Size of each vector /// Array of random vectors - private static Vector[] CreateVectorArray(int arraySize, int vectorSize) + static Vector[] CreateVectorArray(int arraySize, int vectorSize) { var result = new Vector[arraySize]; for (var i = 0; i < result.Length; i++) @@ -710,7 +705,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// Residual data. /// x data. /// b data. - private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { // -Ax = residual matrix.Multiply(x, residual); @@ -728,7 +723,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative /// Source . /// Residual . /// true if continue, otherwise false - private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) + bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) { if (_hasBeenStopped) { @@ -764,7 +759,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative throw new ArgumentNullException("input"); } - var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount); + var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount); Solve(matrix, input, result); return result; } @@ -800,7 +795,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, (Vector)input.Column(column)); + var solution = Solve(matrix, (Vector) input.Column(column)); foreach (var element in solution.GetIndexedEnumerator()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs index a9f125d8..ab9f6680 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs @@ -39,6 +39,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif /// @@ -282,15 +283,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative var temp1 = new DenseVector(input.Count); var temp2 = new DenseVector(input.Count); - // Initialize - var startNorm = input.L2Norm(); - // Define the scalars Complex alpha = 0; Complex eta = 0; double theta = 0; - var tau = startNorm.Real; + // Initialize + var tau = input.L2Norm().Real; Complex rho = tau*tau; // Calculate the initial values for v @@ -311,7 +310,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative if (IsEven(iterationNumber)) { // sigma = (v, r) - var sigma = v.DotProduct(r.Conjugate()); + var sigma = r.ConjugateDotProduct(v); if (sigma.Real.AlmostEqual(0, 1) && sigma.Imaginary.AlmostEqual(0, 1)) { // FAIL HERE @@ -349,7 +348,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative yinternal.Add(temp, d); // theta = ||pseudoResiduals||_2 / tau - theta = pseudoResiduals.L2Norm().Real / tau; + theta = pseudoResiduals.L2Norm().Real/tau; var c = 1/Math.Sqrt(1 + (theta*theta)); // tau = tau * theta * c @@ -392,7 +391,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative break; } - var rhoNew = pseudoResiduals.DotProduct(r.Conjugate()); + var rhoNew = r.ConjugateDotProduct(pseudoResiduals); var beta = rhoNew/rho; // Update rho for the next loop diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs index b23b4b24..fd26c8ad 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs @@ -71,24 +71,24 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// The status used if there is no status, i.e. the solver hasn't run yet and there is no /// iterator. /// - private static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); + static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); /// /// The preconditioner that will be used. Can be set to , in which case the default /// pre-conditioner will be used. /// - private IPreConditioner _preconditioner; + IPreConditioner _preconditioner; /// /// The iterative process controller. /// - private IIterator _iterator; + IIterator _iterator; /// /// Indicates if the user has stopped the solver. /// - private bool _hasBeenStopped; - + bool _hasBeenStopped; + /// /// Initializes a new instance of the class. /// @@ -96,7 +96,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// When using this constructor the solver will use the with /// the standard settings and a default preconditioner. /// - public BiCgStab() : this(null, null) + public BiCgStab() + : this(null, null) { } @@ -119,7 +120,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// /// /// The that will be used to monitor the iterative process. - public BiCgStab(IIterator iterator) : this(null, iterator) + public BiCgStab(IIterator iterator) + : this(null, iterator) { } @@ -131,7 +133,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// the standard settings. /// /// The that will be used to precondition the matrix equation. - public BiCgStab(IPreConditioner preconditioner) : this(preconditioner, null) + public BiCgStab(IPreConditioner preconditioner) + : this(preconditioner, null) { } @@ -181,10 +184,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// public ICalculationStatus IterationResult { - get - { - return (_iterator != null) ? _iterator.Status : DefaultStatus; - } + get { return (_iterator != null) ? _iterator.Status : DefaultStatus; } } /// @@ -273,9 +273,9 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative { _preconditioner = new UnitPreconditioner(); } - + _preconditioner.Initialize(matrix); - + // Compute r_0 = b - Ax_0 for some initial guess x_0 // In this case we take x_0 = vector // This is basically a SAXPY so it could be made a lot faster @@ -307,7 +307,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative { // rho_(i-1) = r~^T r_(i-1) // dotproduct r~ and r_(i-1) var oldRho = currentRho; - currentRho = tempResiduals.DotProduct(residuals); + currentRho = tempResiduals.ConjugateDotProduct(residuals); // if (rho_(i-1) == 0) // METHOD FAILS // If rho is only 1 ULP from zero then we fail. @@ -320,7 +320,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative if (iterationNumber != 0) { // beta_(i-1) = (rho_(i-1)/rho_(i-2))(alpha_(i-1)/omega(i-1)) - var beta = (currentRho / oldRho) * (alpha / omega); + 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, temp); @@ -339,12 +339,12 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative // SOLVE Mp~ = p_i // M = preconditioner _preconditioner.Approximate(vecP, vecPdash); - + // nu_i = Ap~ matrix.Multiply(vecPdash, nu); // alpha_i = rho_(i-1)/ (r~^T nu_i) = rho / dotproduct(r~ and nu_i) - alpha = currentRho * 1 / tempResiduals.DotProduct(nu); + alpha = currentRho*1/tempResiduals.ConjugateDotProduct(nu); // s = r_(i-1) - alpha_i nu_i nu.Multiply(-alpha, temp); @@ -388,7 +388,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative matrix.Multiply(vecSdash, temp); // omega_i = temp^T s / temp^T temp - omega = temp.DotProduct(vecS) / temp.DotProduct(temp); + omega = temp.ConjugateDotProduct(vecS)/temp.ConjugateDotProduct(temp); // x_i = x_(i-1) + alpha_i p^ + omega_i s^ temp.Multiply(-omega, residuals); @@ -432,11 +432,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// Residual values in . /// Instance of the x. /// Instance of the b. - private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { // -Ax = residual matrix.Multiply(x, residual); - + // Do not use residual = residual.Negate() because it creates another object residual.Multiply(-1, residual); @@ -452,7 +452,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// Source . /// Residual . /// true if continue, otherwise false - private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) + bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) { if (_hasBeenStopped) { @@ -488,7 +488,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative throw new ArgumentNullException("input"); } - var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount); + var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount); Solve(matrix, input, result); return result; } @@ -524,7 +524,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, (Vector)input.Column(column)); + var solution = Solve(matrix, (Vector) input.Column(column)); foreach (var element in solution.GetIndexedEnumerator()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs index 26e5fa19..2dd4c991 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs @@ -377,7 +377,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative matrix.Multiply(temp, s); // alpha_k = (r*_0 * r_k) / (r*_0 * s_k) - var alpha = rdash.DotProduct(residuals)/rdash.DotProduct(s); + var alpha = rdash.ConjugateDotProduct(residuals)/rdash.ConjugateDotProduct(s); // y_k = t_(k-1) - r_k - alpha_k * w_(k-1) + alpha_k s_k s.Subtract(w, temp); @@ -399,7 +399,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative // c_k = A d_k matrix.Multiply(temp, c); - var cdot = c.DotProduct(c); + var cdot = c.ConjugateDotProduct(c); // cDot can only be zero if c is a zero vector // We'll set cDot to 1 if it is zero to prevent NaN's @@ -414,7 +414,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative // to do at least one at the start to initialize the // system, but we'll only have to take special measures // if we don't do any so ... - var ctdot = c.DotProduct(t); + var ctdot = c.ConjugateDotProduct(t); Complex32 eta; Complex32 sigma; if (((_numberOfBiCgStabSteps == 0) && (iterationNumber == 0)) || ShouldRunBiCgStabSteps(iterationNumber)) @@ -427,7 +427,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative } else { - var ydot = y.DotProduct(y); + var ydot = y.ConjugateDotProduct(y); // yDot can only be zero if y is a zero vector // We'll set yDot to 1 if it is zero to prevent NaN's @@ -438,8 +438,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative ydot = 1.0f; } - var ytdot = y.DotProduct(t); - var cydot = c.DotProduct(y); + var ytdot = y.ConjugateDotProduct(t); + var cydot = c.ConjugateDotProduct(y); var denom = (cdot*ydot) - (cydot*cydot); @@ -493,7 +493,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative // 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; + beta = (!sigma.Real.AlmostEqual(0, 1) || !sigma.Imaginary.AlmostEqual(0, 1)) ? alpha/sigma*rdash.ConjugateDotProduct(residuals)/rdash.ConjugateDotProduct(t0) : 0; // w_k = c_k + beta_k s_k s.Multiply(beta, temp2); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs index e164fbb9..88d1bd41 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs @@ -68,39 +68,39 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// /// The default number of starting vectors. /// - private const int DefaultNumberOfStartingVectors = 50; - + const int DefaultNumberOfStartingVectors = 50; + /// /// The status used if there is no status, i.e. the solver hasn't run yet and there is no /// iterator. /// - private static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); + static readonly ICalculationStatus DefaultStatus = new CalculationIndetermined(); /// /// The preconditioner that will be used. Can be set to , in which case the default /// pre-conditioner will be used. /// - private IPreConditioner _preconditioner; + IPreConditioner _preconditioner; /// /// The iterative process controller. /// - private IIterator _iterator; + IIterator _iterator; /// /// The collection of starting vectors which are used as the basis for the Krylov sub-space. /// - private IList _startingVectors; + IList _startingVectors; /// /// The number of starting vectors used by the algorithm /// - private int _numberOfStartingVectors = DefaultNumberOfStartingVectors; + int _numberOfStartingVectors = DefaultNumberOfStartingVectors; /// /// Indicates if the user has stopped the solver. /// - private bool _hasBeenStopped; + bool _hasBeenStopped; /// /// Initializes a new instance of the class. @@ -109,7 +109,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// When using this constructor the solver will use the with /// the standard settings and a default preconditioner. /// - public MlkBiCgStab() : this(null, null) + public MlkBiCgStab() + : this(null, null) { } @@ -132,7 +133,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// /// /// The that will be used to monitor the iterative process. - public MlkBiCgStab(IIterator iterator) : this(null, iterator) + public MlkBiCgStab(IIterator iterator) + : this(null, iterator) { } @@ -144,7 +146,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// the standard settings. /// /// The that will be used to precondition the matrix equation. - public MlkBiCgStab(IPreConditioner preconditioner) : this(preconditioner, null) + public MlkBiCgStab(IPreConditioner preconditioner) + : this(preconditioner, null) { } @@ -181,10 +184,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative public int NumberOfStartingVectors { [DebuggerStepThrough] - get - { - return _numberOfStartingVectors; - } + get { return _numberOfStartingVectors; } [DebuggerStepThrough] set @@ -231,10 +231,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative public IList StartingVectors { [DebuggerStepThrough] - get - { - return _startingVectors; - } + get { return _startingVectors; } [DebuggerStepThrough] set @@ -256,10 +253,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative public ICalculationStatus IterationResult { [DebuggerStepThrough] - get - { - return (_iterator != null) ? _iterator.Status : DefaultStatus; - } + get { return (_iterator != null) ? _iterator.Status : DefaultStatus; } } /// @@ -348,7 +342,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative { _preconditioner = new UnitPreconditioner(); } - + _preconditioner.Initialize(matrix); // Choose an initial guess x_0 @@ -402,7 +396,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative Vector zw = new DenseVector(residuals.Count); var d = CreateVectorArray(_startingVectors.Count, residuals.Count); - + // g_0 = r_0 var g = CreateVectorArray(_startingVectors.Count, residuals.Count); residuals.CopyTo(g[k - 1]); @@ -420,14 +414,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative matrix.Multiply(gtemp, w[k - 1]); // c_((j-1)k+k) = q^T_1 w_((j-1)k+k) - c[k - 1] = _startingVectors[0].DotProduct(w[k - 1]); + c[k - 1] = _startingVectors[0].ConjugateDotProduct(w[k - 1]); if (c[k - 1].Real.AlmostEqual(0, 1) && c[k - 1].Imaginary.AlmostEqual(0, 1)) { throw new Exception("Iterative solver experience a numerical break down"); } // alpha_(jk+1) = q^T_1 r_((j-1)k+k) / c_((j-1)k+k) - var alpha = _startingVectors[0].DotProduct(residuals) / c[k - 1]; + var alpha = _startingVectors[0].ConjugateDotProduct(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, temp); @@ -439,7 +433,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative // rho_(j+1) = -u^t_(jk+1) A u~_(jk+1) / ||A u~_(jk+1)||^2 matrix.Multiply(temp1, temp); - var rho = temp.DotProduct(temp); + var rho = temp.ConjugateDotProduct(temp); // If rho is zero then temp is a zero vector and we're probably // about to have zero residuals (i.e. an exact solution). @@ -449,7 +443,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative rho = 1.0f; } - rho = -u.DotProduct(temp) / rho; + rho = -u.ConjugateDotProduct(temp)/rho; // r_(jk+1) = rho_(j+1) A u~_(jk+1) + u_(jk+1) u.CopyTo(residuals); @@ -502,7 +496,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative for (var s = i; s < k - 1; s++) { // beta^(jk+i)_((j-1)k+s) = -q^t_(s+1) z_d / c_((j-1)k+s) - beta = -_startingVectors[s + 1].DotProduct(zd) / c[s]; + beta = -_startingVectors[s + 1].ConjugateDotProduct(zd)/c[s]; // z_d = z_d + beta^(jk+i)_((j-1)k+s) d_((j-1)k+s) d[s].Multiply(beta, temp); @@ -521,7 +515,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative } } - beta = rho * c[k - 1]; + beta = rho*c[k - 1]; if (beta.Real.AlmostEqual(0, 1) && beta.Imaginary.AlmostEqual(0, 1)) { throw new Exception("Iterative solver experience a numerical break down"); @@ -530,7 +524,7 @@ 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)) zw.Multiply(rho, temp2); residuals.Add(temp2, temp); - beta = -_startingVectors[0].DotProduct(temp) / beta; + beta = -_startingVectors[0].ConjugateDotProduct(temp)/beta; // z_g = z_g + beta^(jk+i)_((j-1)k+k) g_((j-1)k+k) g[k - 1].Multiply(beta, temp); @@ -550,7 +544,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative for (var s = 0; s < i - 1; s++) { // beta^(jk+i)_(jk+s) = -q^T_s+1 z_d / c_(jk+s) - beta = -_startingVectors[s + 1].DotProduct(zd) / c[s]; + beta = -_startingVectors[s + 1].ConjugateDotProduct(zd)/c[s]; // z_d = z_d + beta^(jk+i)_(jk+s) * d_(jk+s) d[s].Multiply(beta, temp); @@ -573,14 +567,14 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative if (i < k - 1) { // c_(jk+1) = q^T_i+1 d_(jk+i) - c[i] = _startingVectors[i + 1].DotProduct(d[i]); + c[i] = _startingVectors[i + 1].ConjugateDotProduct(d[i]); if (c[i].Real.AlmostEqual(0, 1) && c[i].Imaginary.AlmostEqual(0, 1)) { throw new Exception("Iterative solver experience a numerical break down"); } // alpha_(jk+i+1) = q^T_(i+1) u_(jk+i) / c_(jk+i) - alpha = _startingVectors[i + 1].DotProduct(u) / c[i]; + alpha = _startingVectors[i + 1].ConjugateDotProduct(u)/c[i]; // u_(jk+i+1) = u_(jk+i) - alpha_(jk+i+1) d_(jk+i) d[i].Multiply(-alpha, temp); @@ -591,7 +585,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative _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, temp); + gtemp.Multiply(rho*alpha, temp); xtemp.Add(temp, temp2); temp2.CopyTo(xtemp); @@ -599,7 +593,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative 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, temp); + w[i].Multiply(-rho*alpha, temp); residuals.Add(temp, temp2); temp2.CopyTo(residuals); @@ -626,7 +620,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// Maximum number /// Number of variables /// Number of starting vectors to create - private static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables) + static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables) { // Create no more starting vectors than the size of the problem - 1 return Math.Min(maximumNumberOfStartingVectors, (numberOfVariables - 1)); @@ -643,7 +637,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// the is smaller than /// the . /// - private static IList CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) + static IList CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) { // Create no more starting vectors than the size of the problem - 1 // Get random values and then orthogonalize them with @@ -662,7 +656,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative var samplesIm = distribution.Samples().Take(matrix.RowCount).ToArray(); for (int j = 0; j < matrix.RowCount; j++) { - samples[j] = new Complex32((float)samplesRe[j], (float)samplesIm[j]); + samples[j] = new Complex32((float) samplesRe[j], (float) samplesIm[j]); } // Set the column @@ -677,10 +671,10 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative var result = new List(); for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) { - result.Add((Vector)orthogonalMatrix.Column(i)); - + result.Add((Vector) orthogonalMatrix.Column(i)); + // Normalize the result vector - result[i].Multiply(1 / result[i].L2Norm().Real, result[i]); + result[i].Multiply(1/result[i].L2Norm().Real, result[i]); } return result; @@ -692,7 +686,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// Number of vectors /// Size of each vector /// Array of random vectors - private static Vector[] CreateVectorArray(int arraySize, int vectorSize) + static Vector[] CreateVectorArray(int arraySize, int vectorSize) { var result = new Vector[arraySize]; for (var i = 0; i < result.Length; i++) @@ -710,7 +704,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// Residual data. /// x data. /// b data. - private static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { // -Ax = residual matrix.Multiply(x, residual); @@ -728,7 +722,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative /// Source . /// Residual . /// true if continue, otherwise false - private bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) + bool ShouldContinue(int iterationNumber, Vector result, Vector source, Vector residuals) { if (_hasBeenStopped) { @@ -764,7 +758,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative throw new ArgumentNullException("input"); } - var result = (Matrix)matrix.CreateMatrix(input.RowCount, input.ColumnCount); + var result = (Matrix) matrix.CreateMatrix(input.RowCount, input.ColumnCount); Solve(matrix, input, result); return result; } @@ -800,7 +794,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative for (var column = 0; column < input.ColumnCount; column++) { - var solution = Solve(matrix, (Vector)input.Column(column)); + var solution = Solve(matrix, (Vector) input.Column(column)); foreach (var element in solution.GetIndexedEnumerator()) { result.At(element.Item1, column, element.Item2); diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs index 64911e68..50619e51 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs @@ -282,15 +282,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative var temp1 = new DenseVector(input.Count); var temp2 = new DenseVector(input.Count); - // Initialize - var startNorm = input.L2Norm(); - // Define the scalars Complex32 alpha = 0; Complex32 eta = 0; float theta = 0; - var tau = startNorm.Real; + // Initialize + var tau = input.L2Norm().Real; Complex32 rho = tau*tau; // Calculate the initial values for v @@ -311,7 +309,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative if (IsEven(iterationNumber)) { // sigma = (v, r) - var sigma = v.DotProduct(r.Conjugate()); + var sigma = r.ConjugateDotProduct(v); if (sigma.Real.AlmostEqual(0, 1) && sigma.Imaginary.AlmostEqual(0, 1)) { // FAIL HERE @@ -349,7 +347,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative yinternal.Add(temp, d); // theta = ||pseudoResiduals||_2 / tau - theta = pseudoResiduals.L2Norm().Real / tau; + theta = pseudoResiduals.L2Norm().Real/tau; var c = 1/(float) Math.Sqrt(1 + (theta*theta)); // tau = tau * theta * c @@ -392,7 +390,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative break; } - var rhoNew = pseudoResiduals.DotProduct(r.Conjugate()); + var rhoNew = r.ConjugateDotProduct(pseudoResiduals); var beta = rhoNew/rho; // Update rho for the next loop diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs index b7b385cf..5177a8e4 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/Iterative/TFQMR.cs @@ -280,17 +280,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers.Iterative var temp = new DenseVector(input.Count); var temp1 = new DenseVector(input.Count); var temp2 = new DenseVector(input.Count); - - // Initialize - var startNorm = input.L2Norm(); // Define the scalars double alpha = 0; double eta = 0; double theta = 0; - var tau = startNorm; - var rho = tau * tau; + // Initialize + var tau = input.L2Norm(); + var rho = tau*tau; // Calculate the initial values for v // M temp = yEven diff --git a/src/Numerics/LinearAlgebra/Generic/Vector.Arithmetic.cs b/src/Numerics/LinearAlgebra/Generic/Vector.Arithmetic.cs index ae74970c..9c475d50 100644 --- a/src/Numerics/LinearAlgebra/Generic/Vector.Arithmetic.cs +++ b/src/Numerics/LinearAlgebra/Generic/Vector.Arithmetic.cs @@ -513,6 +513,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// The sum of a[i]*b[i] for all i. /// If is not of the same size. /// If is . + /// public T DotProduct(Vector other) { if (other == null) throw new ArgumentNullException("other"); @@ -528,6 +529,7 @@ namespace MathNet.Numerics.LinearAlgebra.Generic /// The sum of conj(a[i])*b[i] for all i. /// If is not of the same size. /// If is . + /// public T ConjugateDotProduct(Vector other) { if (other == null) throw new ArgumentNullException("other"); diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs index d1828819..0f6230e2 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs @@ -281,15 +281,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers.Iterative var temp1 = new DenseVector(input.Count); var temp2 = new DenseVector(input.Count); - // Initialize - var startNorm = input.L2Norm(); - // Define the scalars float alpha = 0; float eta = 0; float theta = 0; - var tau = startNorm; + // Initialize + var tau = input.L2Norm(); var rho = tau*tau; // Calculate the initial values for v