diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs index 00faebe7..b32802c2 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/BiCgStab.cs @@ -40,6 +40,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif /// @@ -77,24 +78,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. /// @@ -102,7 +103,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) { } @@ -125,7 +127,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) { } @@ -137,7 +140,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) { } @@ -187,10 +191,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; } } /// @@ -279,9 +280,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 @@ -313,7 +314,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. @@ -326,7 +327,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); @@ -345,12 +346,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); @@ -394,7 +395,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); @@ -438,11 +439,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); @@ -458,7 +459,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) { @@ -494,7 +495,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; } @@ -530,7 +531,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 b1fea710..f5f37b9c 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/GpBiCg.cs @@ -40,6 +40,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif /// @@ -378,7 +379,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); @@ -400,7 +401,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 @@ -415,7 +416,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)) @@ -428,7 +429,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 @@ -439,8 +440,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); @@ -494,7 +495,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 bec572d4..8bf1155d 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/MlkBiCgStab.cs @@ -44,6 +44,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif /// @@ -74,39 +75,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. @@ -115,7 +116,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) { } @@ -138,7 +140,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) { } @@ -150,7 +153,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) { } @@ -187,10 +191,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative public int NumberOfStartingVectors { [DebuggerStepThrough] - get - { - return _numberOfStartingVectors; - } + get { return _numberOfStartingVectors; } [DebuggerStepThrough] set @@ -237,10 +238,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative public IList StartingVectors { [DebuggerStepThrough] - get - { - return _startingVectors; - } + get { return _startingVectors; } [DebuggerStepThrough] set @@ -262,10 +260,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; } } /// @@ -349,7 +344,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative { _preconditioner = new UnitPreconditioner(); } - + _preconditioner.Initialize(matrix); // Choose an initial guess x_0 @@ -403,7 +398,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]); @@ -421,14 +416,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); @@ -440,7 +435,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). @@ -450,7 +445,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); @@ -503,7 +498,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); @@ -522,7 +517,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"); @@ -531,7 +526,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); @@ -551,7 +546,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); @@ -574,14 +569,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); @@ -592,7 +587,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); @@ -600,7 +595,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); @@ -627,7 +622,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)); @@ -644,7 +639,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 @@ -678,10 +673,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; @@ -693,7 +688,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++) @@ -711,7 +706,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); @@ -729,7 +724,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) { @@ -765,7 +760,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; } @@ -801,7 +796,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 1b2f9953..73c7f0b4 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/Iterative/TFQMR.cs @@ -40,6 +40,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers.Iterative using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; + #endif /// @@ -283,15 +284,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 @@ -312,7 +311,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 @@ -350,7 +349,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 @@ -393,7 +392,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 02b3338c..5c738a78 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/BiCgStab.cs @@ -72,24 +72,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. /// @@ -97,7 +97,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) { } @@ -120,7 +121,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) { } @@ -132,7 +134,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) { } @@ -182,10 +185,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; } } /// @@ -274,9 +274,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 @@ -308,7 +308,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. @@ -321,7 +321,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); @@ -340,12 +340,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); @@ -389,7 +389,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); @@ -433,11 +433,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); @@ -453,7 +453,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) { @@ -489,7 +489,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; } @@ -525,7 +525,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 c06af12b..cf19778a 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/GpBiCg.cs @@ -378,7 +378,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); @@ -400,7 +400,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 @@ -415,7 +415,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)) @@ -428,7 +428,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 @@ -439,8 +439,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); @@ -494,7 +494,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 45512f67..70a40c88 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/MlkBiCgStab.cs @@ -69,39 +69,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. @@ -110,7 +110,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) { } @@ -133,7 +134,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) { } @@ -145,7 +147,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) { } @@ -182,10 +185,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative public int NumberOfStartingVectors { [DebuggerStepThrough] - get - { - return _numberOfStartingVectors; - } + get { return _numberOfStartingVectors; } [DebuggerStepThrough] set @@ -232,10 +232,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative public IList StartingVectors { [DebuggerStepThrough] - get - { - return _startingVectors; - } + get { return _startingVectors; } [DebuggerStepThrough] set @@ -257,10 +254,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; } } /// @@ -349,7 +343,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Iterative { _preconditioner = new UnitPreconditioner(); } - + _preconditioner.Initialize(matrix); // Choose an initial guess x_0 @@ -403,7 +397,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]); @@ -421,14 +415,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); @@ -440,7 +434,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). @@ -450,7 +444,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); @@ -503,7 +497,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); @@ -522,7 +516,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"); @@ -531,7 +525,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); @@ -551,7 +545,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); @@ -574,14 +568,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); @@ -592,7 +586,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); @@ -600,7 +594,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); @@ -627,7 +621,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)); @@ -644,7 +638,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 @@ -663,7 +657,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 @@ -678,10 +672,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; @@ -693,7 +687,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++) @@ -711,7 +705,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); @@ -729,7 +723,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) { @@ -765,7 +759,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; } @@ -801,7 +795,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 93b1139a..ef3b8668 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/Iterative/TFQMR.cs @@ -283,15 +283,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 @@ -312,7 +310,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 @@ -350,7 +348,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 @@ -393,7 +391,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 6ca31885..dd8c7421 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/Single/Solvers/Iterative/TFQMR.cs b/src/Numerics/LinearAlgebra/Single/Solvers/Iterative/TFQMR.cs index fdd8c3f0..5d37f29c 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 diff --git a/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs b/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs index a9b37d04..876259ff 100644 --- a/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs +++ b/src/Numerics/LinearAlgebra/Vector.Arithmetic.cs @@ -439,6 +439,7 @@ namespace MathNet.Numerics.LinearAlgebra /// The other vector. /// The sum of a[i]*b[i] for all i. /// If is not of the same size. + /// public T DotProduct(Vector other) { if (Count != other.Count) throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other"); @@ -453,6 +454,7 @@ namespace MathNet.Numerics.LinearAlgebra /// 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 (Count != other.Count) throw new ArgumentException(Resources.ArgumentVectorsSameLength, "other");