diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/BiCgStab.cs index 6cbdc370..76a9da61 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/BiCgStab.cs @@ -73,44 +73,36 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers public sealed class BiCgStab : IIterativeSolver { /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public BiCgStab() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { - } + // -Ax = residual + matrix.Multiply(x, residual); - /// - /// Stops the solve process. - /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() - { - _hasBeenStopped = true; + // Do not use residual = residual.Negate() because it creates another object + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Determine if calculation should continue /// - /// The coefficient , A. - /// The solution , b. - /// The result , x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; } /// @@ -122,32 +114,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The result , x. public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Parameters checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -321,63 +292,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the 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); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient , A. - /// The solution , B. - /// The result , X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -411,5 +325,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient , A. + /// The solution , b. + /// The result , x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient , A. + /// The solution , B. + /// The result , X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs index 55d1f34c..58d72645 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/CompositeSolver.cs @@ -65,57 +65,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// readonly List, IPreconditioner>> _solvers; - /// - /// The status of the calculation. - /// - IterationStatus _status = IterationStatus.Indetermined; - - /// - /// A flag indicating if the solver has been stopped or not. - /// - bool _hasBeenStopped; - - /// - /// The solver that is currently running. Reference is used to be able to stop the - /// solver if the user cancels the solve process. - /// - IIterativeSolver _currentSolver; - public CompositeSolver(IEnumerable> solvers) { _solvers = solvers.Select(setup => new Tuple, IPreconditioner>(setup.CreateSolver(), setup.CreatePreconditioner() ?? new UnitPreconditioner())).ToList(); } - /// - /// Stops the solve process. - /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() - { - _hasBeenStopped = true; - var currentSolver = _currentSolver; - if (currentSolver != null) - { - currentSolver.StopSolve(); - } - } - - /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. - /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. @@ -125,12 +79,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - _currentSolver = null; - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); @@ -153,11 +101,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers var internalInput = input.Clone(); var internalResult = result.Clone(); - foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped)) + foreach (var solver in _solvers) { // Store a reference to the solver so we can stop it. - _currentSolver = solver.Item1; + IterationStatus status; try { // Reset the iterator and pass it to the solver @@ -165,7 +113,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers // Start the solver solver.Item1.Solve(matrix, internalInput, internalResult, iterator, solver.Item2); - _status = iterator.Status; + status = iterator.Status; } catch (Exception) { @@ -174,12 +122,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers // Switch to the next preconditioner. // Reset the solution vector to the previous solution input.CopyTo(internalInput); - _status = IterationStatus.Running; continue; } // There was no fatal breakdown so check the status - if (_status == IterationStatus.Converged) + if (status == IterationStatus.Converged) { // We're done internalResult.CopyTo(result); @@ -189,7 +136,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers // We're not done // Either: // - calculation finished without convergence - if (_status == IterationStatus.StoppedWithoutConvergence) + if (status == IterationStatus.StoppedWithoutConvergence) { // Copy the internal result to the result vector and // continue with the calculation. @@ -203,27 +150,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers input.CopyTo(internalInput); } } - - // Inside the loop we already copied the final results (if there are any) - // So no need to do that again. - - // Clean up - // No longer need the current solver - _currentSolver = null; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; } /// @@ -259,5 +185,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs index 27200050..fb94c6c2 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/GpBiCg.cs @@ -82,22 +82,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// int _numberOfGpbiCgSteps = 4; - /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. - /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public GpBiCg() - { - } - /// /// Gets or sets the number of steps taken with the BiCgStab algorithm /// before switching over to the GPBiCG algorithm. @@ -137,29 +121,51 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } /// - /// Stops the solve process. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually - /// stop the process. - /// - public void StopSolve() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { - _hasBeenStopped = true; + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Determine if calculation should continue /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; + } + + /// + /// Decide if to do steps with BiCgStab + /// + /// Number of iteration + /// true if yes, otherwise false + bool ShouldRunBiCgStabSteps(int iterationNumber) + { + // Run the first steps as BiCGStab + // The number of steps past a whole iteration set + var difference = iterationNumber % (_numberOfBiCgStabSteps + _numberOfGpbiCgSteps); + + // Do steps with BiCGStab if: + // - The difference is zero or more (i.e. we have done zero or more complete cycles) + // - The difference is less than the number of BiCGStab steps that should be taken + return (difference >= 0) && (difference < _numberOfBiCgStabSteps); } /// @@ -171,32 +177,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (input.Count != matrix.RowCount || result.Count != input.Count) { throw Matrix.DimensionsDontMatch(matrix, input, result); @@ -406,78 +391,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the b. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Decide if to do steps with BiCgStab - /// - /// Number of iteration - /// true if yes, otherwise false - bool ShouldRunBiCgStabSteps(int iterationNumber) - { - // Run the first steps as BiCGStab - // The number of steps past a whole iteration set - var difference = iterationNumber%(_numberOfBiCgStabSteps + _numberOfGpbiCgSteps); - - // Do steps with BiCGStab if: - // - The difference is zero or more (i.e. we have done zero or more complete cycles) - // - The difference is less than the number of BiCGStab steps that should be taken - return (difference >= 0) && (difference < _numberOfBiCgStabSteps); - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -511,5 +424,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs index 767de5f8..1bfad10a 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/MlkBiCgStab.cs @@ -32,7 +32,6 @@ using System; using System.Collections.Generic; using System.Diagnostics; using System.Linq; -using System.Threading.Tasks; using MathNet.Numerics.Distributions; using MathNet.Numerics.LinearAlgebra.Solvers; using MathNet.Numerics.Properties; @@ -44,7 +43,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers using Complex = Numerics.Complex; #else using Complex = System.Numerics.Complex; - #endif /// @@ -87,22 +85,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// int _numberOfStartingVectors = DefaultNumberOfStartingVectors; - /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. - /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public MlkBiCgStab() - { - } - /// /// Gets or sets the number of starting vectors. /// @@ -159,30 +141,119 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } /// - /// Stops the solve process. + /// Gets the number of starting vectors to create /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() + /// Maximum number + /// Number of variables + /// Number of starting vectors to create + static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables) { - _hasBeenStopped = true; + // Create no more starting vectors than the size of the problem - 1 + return Math.Min(maximumNumberOfStartingVectors, (numberOfVariables - 1)); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Returns an array of starting vectors. /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// The maximum number of starting vectors that should be created. + /// The number of variables. + /// + /// An array with starting vectors. The array will never be larger than the + /// but it may be smaller if + /// the is smaller than + /// the . + /// + static IList> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); + // Create no more starting vectors than the size of the problem - 1 + // Get random values and then orthogonalize them with + // modified Gramm - Schmidt + var count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables); + + // Get a random set of samples based on the standard normal distribution with + // mean = 0 and sd = 1 + var distribution = new Normal(); + + var matrix = new DenseMatrix(numberOfVariables, count); + for (var i = 0; i < matrix.ColumnCount; i++) + { + var samples = new Complex[matrix.RowCount]; + var samplesRe = distribution.Samples().Take(matrix.RowCount).ToArray(); + var samplesIm = distribution.Samples().Take(matrix.RowCount).ToArray(); + for (int j = 0; j < matrix.RowCount; j++) + { + samples[j] = new Complex(samplesRe[j], samplesIm[j]); + } + + // Set the column + matrix.SetColumn(i, samples); + } + + // Compute the orthogonalization. + var gs = matrix.GramSchmidt(); + var orthogonalMatrix = gs.Q; + + // Now transfer this to vectors + var result = new List>(); + for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) + { + result.Add(orthogonalMatrix.Column(i)); + + // Normalize the result vector + result[i].Multiply(1 / result[i].L2Norm(), result[i]); + } + + return result; + } + + /// + /// Create random vecrors array + /// + /// Number of vectors + /// Size of each vector + /// Array of random vectors + static Vector[] CreateVectorArray(int arraySize, int vectorSize) + { + var result = new Vector[arraySize]; + for (var i = 0; i < result.Length; i++) + { + result[i] = new DenseVector(vectorSize); + } + return result; } + /// + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax + /// + /// Source A. + /// Residual data. + /// x data. + /// b data. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) + { + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); + } + + /// + /// Determine if calculation should continue + /// + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) + { + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; + } + /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. @@ -192,32 +263,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (input.Count != matrix.RowCount || result.Count != input.Count) { throw Matrix.DimensionsDontMatch(matrix, input, result); @@ -506,144 +556,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers xtemp.CopyTo(result); } - /// - /// Gets the number of starting vectors to create - /// - /// Maximum number - /// Number of variables - /// Number of starting vectors to create - 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)); - } - - /// - /// Returns an array of starting vectors. - /// - /// The maximum number of starting vectors that should be created. - /// The number of variables. - /// - /// An array with starting vectors. The array will never be larger than the - /// but it may be smaller if - /// the is smaller than - /// the . - /// - 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 - // modified Gramm - Schmidt - var count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables); - - // Get a random set of samples based on the standard normal distribution with - // mean = 0 and sd = 1 - var distribution = new Normal(); - - var matrix = new DenseMatrix(numberOfVariables, count); - for (var i = 0; i < matrix.ColumnCount; i++) - { - var samples = new Complex[matrix.RowCount]; - var samplesRe = distribution.Samples().Take(matrix.RowCount).ToArray(); - var samplesIm = distribution.Samples().Take(matrix.RowCount).ToArray(); - for (int j = 0; j < matrix.RowCount; j++) - { - samples[j] = new Complex(samplesRe[j], samplesIm[j]); - } - - // Set the column - matrix.SetColumn(i, samples); - } - - // Compute the orthogonalization. - var gs = matrix.GramSchmidt(); - var orthogonalMatrix = gs.Q; - - // Now transfer this to vectors - var result = new List>(); - for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) - { - result.Add(orthogonalMatrix.Column(i)); - - // Normalize the result vector - result[i].Multiply(1/result[i].L2Norm(), result[i]); - } - - return result; - } - - /// - /// Create random vecrors array - /// - /// Number of vectors - /// Size of each vector - /// Array of random vectors - static Vector[] CreateVectorArray(int arraySize, int vectorSize) - { - var result = new Vector[arraySize]; - for (var i = 0; i < result.Length; i++) - { - result[i] = new DenseVector(vectorSize); - } - - return result; - } - - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Source A. - /// Residual data. - /// x data. - /// b data. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -677,5 +589,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs index 357ea2cc..55bd3dc2 100644 --- a/src/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex/Solvers/TFQMR.cs @@ -62,44 +62,44 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers public sealed class TFQMR : IIterativeSolver { /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public TFQMR() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Stops the solve process. + /// Determine if calculation should continue /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - _hasBeenStopped = true; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Is even? /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number to check + /// true if even, otherwise false + static bool IsEven(int number) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + return number % 2 == 0; } /// @@ -111,32 +111,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (input.Count != matrix.RowCount || result.Count != input.Count) { throw Matrix.DimensionsDontMatch(matrix, input, result); @@ -313,71 +292,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the b. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Is even? - /// - /// Number to check - /// true if even, otherwise false - static bool IsEven(int number) - { - return number%2 == 0; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -411,5 +325,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/BiCgStab.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/BiCgStab.cs index a62c6212..8db9969e 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/BiCgStab.cs @@ -67,44 +67,36 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers public sealed class BiCgStab : IIterativeSolver { /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public BiCgStab() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { - } + // -Ax = residual + matrix.Multiply(x, residual); - /// - /// Stops the solve process. - /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() - { - _hasBeenStopped = true; + // Do not use residual = residual.Negate() because it creates another object + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Determine if calculation should continue /// - /// The coefficient , A. - /// The solution , b. - /// The result , x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; } /// @@ -116,32 +108,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The result , x. public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Parameters checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -315,63 +286,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the 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); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient , A. - /// The solution , B. - /// The result , X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -405,5 +319,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient , A. + /// The solution , b. + /// The result , x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient , A. + /// The solution , B. + /// The result , X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/CompositeSolver.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/CompositeSolver.cs index 3f8eb369..75fde196 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/CompositeSolver.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/CompositeSolver.cs @@ -58,57 +58,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// readonly List, IPreconditioner>> _solvers; - /// - /// The status of the calculation. - /// - IterationStatus _status = IterationStatus.Indetermined; - - /// - /// A flag indicating if the solver has been stopped or not. - /// - bool _hasBeenStopped; - - /// - /// The solver that is currently running. Reference is used to be able to stop the - /// solver if the user cancels the solve process. - /// - IIterativeSolver _currentSolver; - public CompositeSolver(IEnumerable> solvers) { _solvers = solvers.Select(setup => new Tuple, IPreconditioner>(setup.CreateSolver(), setup.CreatePreconditioner() ?? new UnitPreconditioner())).ToList(); } - /// - /// Stops the solve process. - /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() - { - _hasBeenStopped = true; - var currentSolver = _currentSolver; - if (currentSolver != null) - { - currentSolver.StopSolve(); - } - } - - /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. - /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. @@ -118,12 +72,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - _currentSolver = null; - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); @@ -146,11 +94,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers var internalInput = input.Clone(); var internalResult = result.Clone(); - foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped)) + foreach (var solver in _solvers) { // Store a reference to the solver so we can stop it. - _currentSolver = solver.Item1; + IterationStatus status; try { // Reset the iterator and pass it to the solver @@ -158,7 +106,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers // Start the solver solver.Item1.Solve(matrix, internalInput, internalResult, iterator, solver.Item2); - _status = iterator.Status; + status = iterator.Status; } catch (Exception) { @@ -167,12 +115,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers // Switch to the next preconditioner. // Reset the solution vector to the previous solution input.CopyTo(internalInput); - _status = IterationStatus.Running; continue; } // There was no fatal breakdown so check the status - if (_status == IterationStatus.Converged) + if (status == IterationStatus.Converged) { // We're done internalResult.CopyTo(result); @@ -182,7 +129,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers // We're not done // Either: // - calculation finished without convergence - if (_status == IterationStatus.StoppedWithoutConvergence) + if (status == IterationStatus.StoppedWithoutConvergence) { // Copy the internal result to the result vector and // continue with the calculation. @@ -196,27 +143,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers input.CopyTo(internalInput); } } - - // Inside the loop we already copied the final results (if there are any) - // So no need to do that again. - - // Clean up - // No longer need the current solver - _currentSolver = null; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; } /// @@ -252,5 +178,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/GpBiCg.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/GpBiCg.cs index 3c5ae9b2..a04fcce4 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/GpBiCg.cs @@ -75,22 +75,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// int _numberOfGpbiCgSteps = 4; - /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. - /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public GpBiCg() - { - } - /// /// Gets or sets the number of steps taken with the BiCgStab algorithm /// before switching over to the GPBiCG algorithm. @@ -130,29 +114,51 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } /// - /// Stops the solve process. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually - /// stop the process. - /// - public void StopSolve() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { - _hasBeenStopped = true; + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Determine if calculation should continue /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; + } + + /// + /// Decide if to do steps with BiCgStab + /// + /// Number of iteration + /// true if yes, otherwise false + bool ShouldRunBiCgStabSteps(int iterationNumber) + { + // Run the first steps as BiCGStab + // The number of steps past a whole iteration set + var difference = iterationNumber % (_numberOfBiCgStabSteps + _numberOfGpbiCgSteps); + + // Do steps with BiCGStab if: + // - The difference is zero or more (i.e. we have done zero or more complete cycles) + // - The difference is less than the number of BiCGStab steps that should be taken + return (difference >= 0) && (difference < _numberOfBiCgStabSteps); } /// @@ -164,32 +170,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -404,78 +389,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the b. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Decide if to do steps with BiCgStab - /// - /// Number of iteration - /// true if yes, otherwise false - bool ShouldRunBiCgStabSteps(int iterationNumber) - { - // Run the first steps as BiCGStab - // The number of steps past a whole iteration set - var difference = iterationNumber%(_numberOfBiCgStabSteps + _numberOfGpbiCgSteps); - - // Do steps with BiCGStab if: - // - The difference is zero or more (i.e. we have done zero or more complete cycles) - // - The difference is less than the number of BiCGStab steps that should be taken - return (difference >= 0) && (difference < _numberOfBiCgStabSteps); - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -509,5 +422,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/MlkBiCgStab.cs index 4c44c535..4ecea50f 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/MlkBiCgStab.cs @@ -33,7 +33,6 @@ using System.Collections.Generic; using System.Diagnostics; using System.Linq; using MathNet.Numerics.Distributions; -using MathNet.Numerics.LinearAlgebra.Complex32.Solvers.Preconditioners; using MathNet.Numerics.LinearAlgebra.Solvers; using MathNet.Numerics.Properties; @@ -79,22 +78,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// int _numberOfStartingVectors = DefaultNumberOfStartingVectors; - /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. - /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public MlkBiCgStab() - { - } - /// /// Gets or sets the number of starting vectors. /// @@ -151,30 +134,119 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } /// - /// Stops the solve process. + /// Gets the number of starting vectors to create /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() + /// Maximum number + /// Number of variables + /// Number of starting vectors to create + static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables) { - _hasBeenStopped = true; + // Create no more starting vectors than the size of the problem - 1 + return Math.Min(maximumNumberOfStartingVectors, (numberOfVariables - 1)); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Returns an array of starting vectors. /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// The maximum number of starting vectors that should be created. + /// The number of variables. + /// + /// An array with starting vectors. The array will never be larger than the + /// but it may be smaller if + /// the is smaller than + /// the . + /// + static IList> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); + // Create no more starting vectors than the size of the problem - 1 + // Get random values and then orthogonalize them with + // modified Gramm - Schmidt + var count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables); + + // Get a random set of samples based on the standard normal distribution with + // mean = 0 and sd = 1 + var distribution = new Normal(); + + var matrix = new DenseMatrix(numberOfVariables, count); + for (var i = 0; i < matrix.ColumnCount; i++) + { + var samples = new Numerics.Complex32[matrix.RowCount]; + var samplesRe = distribution.Samples().Take(matrix.RowCount).ToArray(); + var samplesIm = distribution.Samples().Take(matrix.RowCount).ToArray(); + for (int j = 0; j < matrix.RowCount; j++) + { + samples[j] = new Numerics.Complex32((float)samplesRe[j], (float)samplesIm[j]); + } + + // Set the column + matrix.SetColumn(i, samples); + } + + // Compute the orthogonalization. + var gs = matrix.GramSchmidt(); + var orthogonalMatrix = gs.Q; + + // Now transfer this to vectors + var result = new List>(); + for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) + { + result.Add(orthogonalMatrix.Column(i)); + + // Normalize the result vector + result[i].Multiply(1 / result[i].L2Norm().Real, result[i]); + } + + return result; + } + + /// + /// Create random vectors array + /// + /// Number of vectors + /// Size of each vector + /// Array of random vectors + static Vector[] CreateVectorArray(int arraySize, int vectorSize) + { + var result = new Vector[arraySize]; + for (var i = 0; i < result.Length; i++) + { + result[i] = new DenseVector(vectorSize); + } + return result; } + /// + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax + /// + /// Source A. + /// Residual data. + /// x data. + /// b data. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) + { + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); + } + + /// + /// Determine if calculation should continue + /// + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) + { + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; + } + /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. @@ -184,32 +256,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -503,144 +554,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers xtemp.CopyTo(result); } - /// - /// Gets the number of starting vectors to create - /// - /// Maximum number - /// Number of variables - /// Number of starting vectors to create - 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)); - } - - /// - /// Returns an array of starting vectors. - /// - /// The maximum number of starting vectors that should be created. - /// The number of variables. - /// - /// An array with starting vectors. The array will never be larger than the - /// but it may be smaller if - /// the is smaller than - /// the . - /// - 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 - // modified Gramm - Schmidt - var count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables); - - // Get a random set of samples based on the standard normal distribution with - // mean = 0 and sd = 1 - var distribution = new Normal(); - - var matrix = new DenseMatrix(numberOfVariables, count); - for (var i = 0; i < matrix.ColumnCount; i++) - { - var samples = new Numerics.Complex32[matrix.RowCount]; - var samplesRe = distribution.Samples().Take(matrix.RowCount).ToArray(); - var samplesIm = distribution.Samples().Take(matrix.RowCount).ToArray(); - for (int j = 0; j < matrix.RowCount; j++) - { - samples[j] = new Numerics.Complex32((float) samplesRe[j], (float) samplesIm[j]); - } - - // Set the column - matrix.SetColumn(i, samples); - } - - // Compute the orthogonalization. - var gs = matrix.GramSchmidt(); - var orthogonalMatrix = gs.Q; - - // Now transfer this to vectors - var result = new List>(); - for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) - { - result.Add(orthogonalMatrix.Column(i)); - - // Normalize the result vector - result[i].Multiply(1/result[i].L2Norm().Real, result[i]); - } - - return result; - } - - /// - /// Create random vectors array - /// - /// Number of vectors - /// Size of each vector - /// Array of random vectors - static Vector[] CreateVectorArray(int arraySize, int vectorSize) - { - var result = new Vector[arraySize]; - for (var i = 0; i < result.Length; i++) - { - result[i] = new DenseVector(vectorSize); - } - - return result; - } - - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Source A. - /// Residual data. - /// x data. - /// b data. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -674,5 +587,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Complex32/Solvers/TFQMR.cs b/src/Numerics/LinearAlgebra/Complex32/Solvers/TFQMR.cs index b11ca1c8..d024a57b 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Solvers/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Solvers/TFQMR.cs @@ -54,44 +54,44 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers public sealed class TFQMR : IIterativeSolver { /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public TFQMR() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Stops the solve process. + /// Determine if calculation should continue /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - _hasBeenStopped = true; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Is even? /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number to check + /// true if even, otherwise false + static bool IsEven(int number) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + return number % 2 == 0; } /// @@ -103,32 +103,11 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -310,71 +289,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the b. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Is even? - /// - /// Number to check - /// true if even, otherwise false - static bool IsEven(int number) - { - return number%2 == 0; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -408,5 +322,33 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/BiCgStab.cs b/src/Numerics/LinearAlgebra/Double/Solvers/BiCgStab.cs index b2f9ad55..306ebff9 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/BiCgStab.cs @@ -66,44 +66,36 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers public sealed class BiCgStab : IIterativeSolver { /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public BiCgStab() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { - } + // -Ax = residual + matrix.Multiply(x, residual); - /// - /// Stops the solve process. - /// - /// - /// It may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() - { - _hasBeenStopped = true; + // Do not use residual = residual.Negate() because it creates another object + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Determine if calculation should continue /// - /// The coefficient , A. - /// The solution , b. - /// The result , x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; } /// @@ -115,32 +107,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The result , x. public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Parameters checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -314,63 +285,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the 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); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient , A. - /// The solution , B. - /// The result , X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -404,5 +318,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient , A. + /// The solution , b. + /// The result , x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient , A. + /// The solution , B. + /// The result , X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/CompositeSolver.cs b/src/Numerics/LinearAlgebra/Double/Solvers/CompositeSolver.cs index 0183ce12..b26fad67 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/CompositeSolver.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/CompositeSolver.cs @@ -58,57 +58,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// readonly List, IPreconditioner>> _solvers; - /// - /// The status of the calculation. - /// - IterationStatus _status = IterationStatus.Indetermined; - - /// - /// A flag indicating if the solver has been stopped or not. - /// - bool _hasBeenStopped; - - /// - /// The solver that is currently running. Reference is used to be able to stop the - /// solver if the user cancels the solve process. - /// - IIterativeSolver _currentSolver; - public CompositeSolver(IEnumerable> solvers) { _solvers = solvers.Select(setup => new Tuple, IPreconditioner>(setup.CreateSolver(), setup.CreatePreconditioner() ?? new UnitPreconditioner())).ToList(); } - /// - /// Stops the solve process. - /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() - { - _hasBeenStopped = true; - var currentSolver = _currentSolver; - if (currentSolver != null) - { - currentSolver.StopSolve(); - } - } - - /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. - /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. @@ -118,12 +72,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - _currentSolver = null; - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); @@ -146,11 +94,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers var internalInput = input.Clone(); var internalResult = result.Clone(); - foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped)) + foreach (var solver in _solvers) { // Store a reference to the solver so we can stop it. - _currentSolver = solver.Item1; + IterationStatus status; try { // Reset the iterator and pass it to the solver @@ -158,7 +106,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers // Start the solver solver.Item1.Solve(matrix, internalInput, internalResult, iterator, solver.Item2); - _status = iterator.Status; + status = iterator.Status; } catch (Exception) { @@ -167,12 +115,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers // Switch to the next preconditioner. // Reset the solution vector to the previous solution input.CopyTo(internalInput); - _status = IterationStatus.Running; continue; } // There was no fatal breakdown so check the status - if (_status == IterationStatus.Converged) + if (status == IterationStatus.Converged) { // We're done internalResult.CopyTo(result); @@ -182,7 +129,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers // We're not done // Either: // - calculation finished without convergence - if (_status == IterationStatus.StoppedWithoutConvergence) + if (status == IterationStatus.StoppedWithoutConvergence) { // Copy the internal result to the result vector and // continue with the calculation. @@ -196,27 +143,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers input.CopyTo(internalInput); } } - - // Inside the loop we already copied the final results (if there are any) - // So no need to do that again. - - // Clean up - // No longer need the current solver - _currentSolver = null; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; } /// @@ -252,5 +178,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/GpBiCg.cs b/src/Numerics/LinearAlgebra/Double/Solvers/GpBiCg.cs index d55612da..634de3c6 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/GpBiCg.cs @@ -75,22 +75,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// int _numberOfGpbiCgSteps = 4; - /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. - /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public GpBiCg() - { - } - /// /// Gets or sets the number of steps taken with the BiCgStab algorithm /// before switching over to the GPBiCG algorithm. @@ -136,29 +120,51 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } /// - /// Stops the solve process. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually - /// stop the process. - /// - public void StopSolve() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { - _hasBeenStopped = true; + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Determine if calculation should continue /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; + } + + /// + /// Decide if to do steps with BiCgStab + /// + /// Number of iteration + /// true if yes, otherwise false + bool ShouldRunBiCgStabSteps(int iterationNumber) + { + // Run the first steps as BiCGStab + // The number of steps past a whole iteration set + var difference = iterationNumber % (_numberOfBiCgStabSteps + _numberOfGpbiCgSteps); + + // Do steps with BiCGStab if: + // - The difference is zero or more (i.e. we have done zero or more complete cycles) + // - The difference is less than the number of BiCGStab steps that should be taken + return (difference >= 0) && (difference < _numberOfBiCgStabSteps); } /// @@ -170,32 +176,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -410,78 +395,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the b. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Decide if to do steps with BiCgStab - /// - /// Number of iteration - /// true if yes, otherwise false - bool ShouldRunBiCgStabSteps(int iterationNumber) - { - // Run the first steps as BiCGStab - // The number of steps past a whole iteration set - var difference = iterationNumber%(_numberOfBiCgStabSteps + _numberOfGpbiCgSteps); - - // Do steps with BiCGStab if: - // - The difference is zero or more (i.e. we have done zero or more complete cycles) - // - The difference is less than the number of BiCGStab steps that should be taken - return (difference >= 0) && (difference < _numberOfBiCgStabSteps); - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -515,5 +428,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Double/Solvers/MlkBiCgStab.cs index 8b650f86..a9fba4ef 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/MlkBiCgStab.cs @@ -78,22 +78,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// int _numberOfStartingVectors = DefaultNumberOfStartingVectors; - /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. - /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public MlkBiCgStab() - { - } - /// /// Gets or sets the number of starting vectors. /// @@ -156,30 +140,113 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } /// - /// Stops the solve process. + /// Gets the number of starting vectors to create /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() + /// Maximum number + /// Number of variables + /// Number of starting vectors to create + static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables) { - _hasBeenStopped = true; + // Create no more starting vectors than the size of the problem - 1 + return Math.Min(maximumNumberOfStartingVectors, (numberOfVariables - 1)); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Returns an array of starting vectors. /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// The maximum number of starting vectors that should be created. + /// The number of variables. + /// + /// An array with starting vectors. The array will never be larger than the + /// but it may be smaller if + /// the is smaller than + /// the . + /// + static IList> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); + // Create no more starting vectors than the size of the problem - 1 + // Get random values and then orthogonalize them with + // modified Gramm - Schmidt + var count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables); + + // Get a random set of samples based on the standard normal distribution with + // mean = 0 and sd = 1 + var distribution = new Normal(); + + var matrix = new DenseMatrix(numberOfVariables, count); + for (var i = 0; i < matrix.ColumnCount; i++) + { + var samples = distribution.Samples().Take(matrix.RowCount).ToArray(); + + // Set the column + matrix.SetColumn(i, samples); + } + + // Compute the orthogonalization. + var gs = matrix.GramSchmidt(); + var orthogonalMatrix = gs.Q; + + // Now transfer this to vectors + var result = new List>(); + for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) + { + result.Add(orthogonalMatrix.Column(i)); + + // Normalize the result vector + result[i].Multiply(1 / result[i].L2Norm(), result[i]); + } + return result; } + /// + /// Create random vecrors array + /// + /// Number of vectors + /// Size of each vector + /// Array of random vectors + static Vector[] CreateVectorArray(int arraySize, int vectorSize) + { + var result = new Vector[arraySize]; + for (var i = 0; i < result.Length; i++) + { + result[i] = new DenseVector(vectorSize); + } + + return result; + } + + /// + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax + /// + /// Source A. + /// Residual data. + /// x data. + /// b data. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) + { + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); + } + + /// + /// Determine if calculation should continue + /// + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) + { + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; + } + /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. @@ -189,32 +256,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -508,138 +554,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers xtemp.CopyTo(result); } - /// - /// Gets the number of starting vectors to create - /// - /// Maximum number - /// Number of variables - /// Number of starting vectors to create - 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)); - } - - /// - /// Returns an array of starting vectors. - /// - /// The maximum number of starting vectors that should be created. - /// The number of variables. - /// - /// An array with starting vectors. The array will never be larger than the - /// but it may be smaller if - /// the is smaller than - /// the . - /// - 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 - // modified Gramm - Schmidt - var count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables); - - // Get a random set of samples based on the standard normal distribution with - // mean = 0 and sd = 1 - var distribution = new Normal(); - - var matrix = new DenseMatrix(numberOfVariables, count); - for (var i = 0; i < matrix.ColumnCount; i++) - { - var samples = distribution.Samples().Take(matrix.RowCount).ToArray(); - - // Set the column - matrix.SetColumn(i, samples); - } - - // Compute the orthogonalization. - var gs = matrix.GramSchmidt(); - var orthogonalMatrix = gs.Q; - - // Now transfer this to vectors - var result = new List>(); - for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) - { - result.Add(orthogonalMatrix.Column(i)); - - // Normalize the result vector - result[i].Multiply(1/result[i].L2Norm(), result[i]); - } - - return result; - } - - /// - /// Create random vecrors array - /// - /// Number of vectors - /// Size of each vector - /// Array of random vectors - static Vector[] CreateVectorArray(int arraySize, int vectorSize) - { - var result = new Vector[arraySize]; - for (var i = 0; i < result.Length; i++) - { - result[i] = new DenseVector(vectorSize); - } - - return result; - } - - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Source A. - /// Residual data. - /// x data. - /// b data. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -673,5 +587,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Double/Solvers/TFQMR.cs b/src/Numerics/LinearAlgebra/Double/Solvers/TFQMR.cs index 40589766..498dcd7f 100644 --- a/src/Numerics/LinearAlgebra/Double/Solvers/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Double/Solvers/TFQMR.cs @@ -54,44 +54,44 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers public sealed class TFQMR : IIterativeSolver { /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public TFQMR() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Stops the solve process. + /// Determine if calculation should continue /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - _hasBeenStopped = true; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Is even? /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number to check + /// true if even, otherwise false + static bool IsEven(int number) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + return number % 2 == 0; } /// @@ -103,32 +103,11 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -310,71 +289,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the b. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Is even? - /// - /// Number to check - /// true if even, otherwise false - static bool IsEven(int number) - { - return number%2 == 0; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -408,5 +322,33 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/BiCgStab.cs b/src/Numerics/LinearAlgebra/Single/Solvers/BiCgStab.cs index 486b69a9..0b48e54f 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/BiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/BiCgStab.cs @@ -66,44 +66,36 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers public sealed class BiCgStab : IIterativeSolver { /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public BiCgStab() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { - } + // -Ax = residual + matrix.Multiply(x, residual); - /// - /// Stops the solve process. - /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() - { - _hasBeenStopped = true; + // Do not use residual = residual.Negate() because it creates another object + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Determine if calculation should continue /// - /// The coefficient , A. - /// The solution , b. - /// The result , x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; } /// @@ -115,32 +107,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The result , x. public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Parameters checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -314,63 +285,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the 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); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient , A. - /// The solution , B. - /// The result , X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -404,5 +318,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient , A. + /// The solution , b. + /// The result , x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient , A. + /// The solution , B. + /// The result , X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/CompositeSolver.cs b/src/Numerics/LinearAlgebra/Single/Solvers/CompositeSolver.cs index 625515d8..c7bee637 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/CompositeSolver.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/CompositeSolver.cs @@ -58,57 +58,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// readonly List, IPreconditioner>> _solvers; - /// - /// The status of the calculation. - /// - IterationStatus _status = IterationStatus.Indetermined; - - /// - /// A flag indicating if the solver has been stopped or not. - /// - bool _hasBeenStopped; - - /// - /// The solver that is currently running. Reference is used to be able to stop the - /// solver if the user cancels the solve process. - /// - IIterativeSolver _currentSolver; - public CompositeSolver(IEnumerable> solvers) { _solvers = solvers.Select(setup => new Tuple, IPreconditioner>(setup.CreateSolver(), setup.CreatePreconditioner() ?? new UnitPreconditioner())).ToList(); } - /// - /// Stops the solve process. - /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() - { - _hasBeenStopped = true; - var currentSolver = _currentSolver; - if (currentSolver != null) - { - currentSolver.StopSolve(); - } - } - - /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. - /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. @@ -118,12 +72,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - _currentSolver = null; - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); @@ -146,11 +94,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers var internalInput = input.Clone(); var internalResult = result.Clone(); - foreach (var solver in _solvers.TakeWhile(solver => !_hasBeenStopped)) + foreach (var solver in _solvers) { // Store a reference to the solver so we can stop it. - _currentSolver = solver.Item1; + IterationStatus status; try { // Reset the iterator and pass it to the solver @@ -158,7 +106,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers // Start the solver solver.Item1.Solve(matrix, internalInput, internalResult, iterator, solver.Item2); - _status = iterator.Status; + status = iterator.Status; } catch (Exception) { @@ -167,12 +115,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers // Switch to the next preconditioner. // Reset the solution vector to the previous solution input.CopyTo(internalInput); - _status = IterationStatus.Running; continue; } // There was no fatal breakdown so check the status - if (_status == IterationStatus.Converged) + if (status == IterationStatus.Converged) { // We're done internalResult.CopyTo(result); @@ -182,7 +129,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers // We're not done // Either: // - calculation finished without convergence - if (_status == IterationStatus.StoppedWithoutConvergence) + if (status == IterationStatus.StoppedWithoutConvergence) { // Copy the internal result to the result vector and // continue with the calculation. @@ -196,27 +143,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers input.CopyTo(internalInput); } } - - // Inside the loop we already copied the final results (if there are any) - // So no need to do that again. - - // Clean up - // No longer need the current solver - _currentSolver = null; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; } /// @@ -252,5 +178,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/GpBiCg.cs b/src/Numerics/LinearAlgebra/Single/Solvers/GpBiCg.cs index 8a8c5046..241c84b2 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/GpBiCg.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/GpBiCg.cs @@ -75,22 +75,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// int _numberOfGpbiCgSteps = 4; - /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. - /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public GpBiCg() - { - } - /// /// Gets or sets the number of steps taken with the BiCgStab algorithm /// before switching over to the GPBiCG algorithm. @@ -130,29 +114,51 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } /// - /// Stops the solve process. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually - /// stop the process. - /// - public void StopSolve() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { - _hasBeenStopped = true; + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Determine if calculation should continue /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; + } + + /// + /// Decide if to do steps with BiCgStab + /// + /// Number of iteration + /// true if yes, otherwise false + bool ShouldRunBiCgStabSteps(int iterationNumber) + { + // Run the first steps as BiCGStab + // The number of steps past a whole iteration set + var difference = iterationNumber % (_numberOfBiCgStabSteps + _numberOfGpbiCgSteps); + + // Do steps with BiCGStab if: + // - The difference is zero or more (i.e. we have done zero or more complete cycles) + // - The difference is less than the number of BiCGStab steps that should be taken + return (difference >= 0) && (difference < _numberOfBiCgStabSteps); } /// @@ -164,32 +170,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -404,78 +389,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the b. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Decide if to do steps with BiCgStab - /// - /// Number of iteration - /// true if yes, otherwise false - bool ShouldRunBiCgStabSteps(int iterationNumber) - { - // Run the first steps as BiCGStab - // The number of steps past a whole iteration set - var difference = iterationNumber%(_numberOfBiCgStabSteps + _numberOfGpbiCgSteps); - - // Do steps with BiCGStab if: - // - The difference is zero or more (i.e. we have done zero or more complete cycles) - // - The difference is less than the number of BiCGStab steps that should be taken - return (difference >= 0) && (difference < _numberOfBiCgStabSteps); - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -509,5 +422,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/MlkBiCgStab.cs b/src/Numerics/LinearAlgebra/Single/Solvers/MlkBiCgStab.cs index 11a5532d..8cc1a5ae 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/MlkBiCgStab.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/MlkBiCgStab.cs @@ -77,22 +77,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// int _numberOfStartingVectors = DefaultNumberOfStartingVectors; - /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. - /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public MlkBiCgStab() - { - } - /// /// Gets or sets the number of starting vectors. /// @@ -155,30 +139,117 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } /// - /// Stops the solve process. + /// Gets the number of starting vectors to create /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() + /// Maximum number + /// Number of variables + /// Number of starting vectors to create + static int NumberOfStartingVectorsToCreate(int maximumNumberOfStartingVectors, int numberOfVariables) { - _hasBeenStopped = true; + // Create no more starting vectors than the size of the problem - 1 + return Math.Min(maximumNumberOfStartingVectors, (numberOfVariables - 1)); } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Returns an array of starting vectors. /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// The maximum number of starting vectors that should be created. + /// The number of variables. + /// + /// An array with starting vectors. The array will never be larger than the + /// but it may be smaller if + /// the is smaller than + /// the . + /// + static IList> CreateStartingVectors(int maximumNumberOfStartingVectors, int numberOfVariables) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); + // Create no more starting vectors than the size of the problem - 1 + // Get random values and then orthogonalize them with + // modified Gramm - Schmidt + var count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables); + + // Get a random set of samples based on the standard normal distribution with + // mean = 0 and sd = 1 + var distribution = new Normal(); + + var matrix = new DenseMatrix(numberOfVariables, count); + for (var i = 0; i < matrix.ColumnCount; i++) + { + var samples = new float[matrix.RowCount]; + for (var j = 0; j < matrix.RowCount; j++) + { + samples[j] = (float)distribution.Sample(); + } + + // Set the column + matrix.SetColumn(i, samples); + } + + // Compute the orthogonalization. + var gs = matrix.GramSchmidt(); + var orthogonalMatrix = gs.Q; + + // Now transfer this to vectors + var result = new List>(); + for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) + { + result.Add(orthogonalMatrix.Column(i)); + + // Normalize the result vector + result[i].Multiply(1 / result[i].L2Norm(), result[i]); + } + return result; } + /// + /// Create random vectors array + /// + /// Number of vectors + /// Size of each vector + /// Array of random vectors + static Vector[] CreateVectorArray(int arraySize, int vectorSize) + { + var result = new Vector[arraySize]; + for (var i = 0; i < result.Length; i++) + { + result[i] = new DenseVector(vectorSize); + } + + return result; + } + + /// + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax + /// + /// Source A. + /// Residual data. + /// x data. + /// b data. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) + { + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); + } + + /// + /// Determine if calculation should continue + /// + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) + { + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; + } + /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector. @@ -188,32 +259,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -507,142 +557,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers xtemp.CopyTo(result); } - /// - /// Gets the number of starting vectors to create - /// - /// Maximum number - /// Number of variables - /// Number of starting vectors to create - 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)); - } - - /// - /// Returns an array of starting vectors. - /// - /// The maximum number of starting vectors that should be created. - /// The number of variables. - /// - /// An array with starting vectors. The array will never be larger than the - /// but it may be smaller if - /// the is smaller than - /// the . - /// - 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 - // modified Gramm - Schmidt - var count = NumberOfStartingVectorsToCreate(maximumNumberOfStartingVectors, numberOfVariables); - - // Get a random set of samples based on the standard normal distribution with - // mean = 0 and sd = 1 - var distribution = new Normal(); - - var matrix = new DenseMatrix(numberOfVariables, count); - for (var i = 0; i < matrix.ColumnCount; i++) - { - var samples = new float[matrix.RowCount]; - for (var j = 0; j < matrix.RowCount; j++) - { - samples[j] = (float) distribution.Sample(); - } - - // Set the column - matrix.SetColumn(i, samples); - } - - // Compute the orthogonalization. - var gs = matrix.GramSchmidt(); - var orthogonalMatrix = gs.Q; - - // Now transfer this to vectors - var result = new List>(); - for (var i = 0; i < orthogonalMatrix.ColumnCount; i++) - { - result.Add(orthogonalMatrix.Column(i)); - - // Normalize the result vector - result[i].Multiply(1/result[i].L2Norm(), result[i]); - } - - return result; - } - - /// - /// Create random vectors array - /// - /// Number of vectors - /// Size of each vector - /// Array of random vectors - static Vector[] CreateVectorArray(int arraySize, int vectorSize) - { - var result = new Vector[arraySize]; - for (var i = 0; i < result.Length; i++) - { - result[i] = new DenseVector(vectorSize); - } - - return result; - } - - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Source A. - /// Residual data. - /// x data. - /// b data. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -676,5 +590,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Single/Solvers/TFQMR.cs b/src/Numerics/LinearAlgebra/Single/Solvers/TFQMR.cs index d3ece71a..26a596fe 100644 --- a/src/Numerics/LinearAlgebra/Single/Solvers/TFQMR.cs +++ b/src/Numerics/LinearAlgebra/Single/Solvers/TFQMR.cs @@ -54,44 +54,44 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers public sealed class TFQMR : IIterativeSolver { /// - /// Indicates if the user has stopped the solver. - /// - bool _hasBeenStopped; - - /// - /// Initializes a new instance of the class. + /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax /// - /// - /// When using this constructor the solver will use the with - /// the standard settings and a default preconditioner. - /// - public TFQMR() + /// Instance of the A. + /// Residual values in . + /// Instance of the x. + /// Instance of the b. + static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) { + // -Ax = residual + matrix.Multiply(x, residual); + residual.Multiply(-1, residual); + + // residual + b + residual.Add(b, residual); } /// - /// Stops the solve process. + /// Determine if calculation should continue /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - public void StopSolve() + /// Number of iterations passed + /// Result . + /// Source . + /// Residual . + /// true if continue, otherwise false + static bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) { - _hasBeenStopped = true; + var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); + return status == IterationStatus.Running || status == IterationStatus.Indetermined; } /// - /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the - /// solution vector and x is the unknown vector. + /// Is even? /// - /// The coefficient matrix, A. - /// The solution vector, b. - /// The result vector, x. - public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + /// Number to check + /// true if even, otherwise false + static bool IsEven(int number) { - var result = new DenseVector(matrix.RowCount); - Solve(matrix, vector, result, iterator, preconditioner); - return result; + return number % 2 == 0; } /// @@ -103,32 +103,11 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers /// The result vector, x public void Solve(Matrix matrix, Vector input, Vector result, Iterator iterator = null, IPreconditioner preconditioner = null) { - // If we were stopped before, we are no longer - // We're doing this at the start of the method to ensure - // that we can use these fields immediately. - _hasBeenStopped = false; - - // Error checks - if (matrix == null) - { - throw new ArgumentNullException("matrix"); - } - if (matrix.RowCount != matrix.ColumnCount) { throw new ArgumentException(Resources.ArgumentMatrixSquare, "matrix"); } - if (input == null) - { - throw new ArgumentNullException("input"); - } - - if (result == null) - { - throw new ArgumentNullException("result"); - } - if (result.Count != input.Count) { throw new ArgumentException(Resources.ArgumentVectorsSameLength); @@ -310,71 +289,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } } - /// - /// Calculates the true residual of the matrix equation Ax = b according to: residual = b - Ax - /// - /// Instance of the A. - /// Residual values in . - /// Instance of the x. - /// Instance of the b. - static void CalculateTrueResidual(Matrix matrix, Vector residual, Vector x, Vector b) - { - // -Ax = residual - matrix.Multiply(x, residual); - residual.Multiply(-1, residual); - - // residual + b - residual.Add(b, residual); - } - - /// - /// Determine if calculation should continue - /// - /// Number of iterations passed - /// Result . - /// Source . - /// Residual . - /// true if continue, otherwise false - bool ShouldContinue(Iterator iterator, int iterationNumber, Vector result, Vector source, Vector residuals) - { - // We stop if either: - // - the user has stopped the calculation - // - the calculation needs to be stopped from a numerical point of view (divergence, convergence etc.) - - if (_hasBeenStopped) - { - iterator.Cancel(); - return true; - } - - var status = iterator.DetermineStatus(iterationNumber, result, source, residuals); - return status == IterationStatus.Running || status == IterationStatus.Indetermined; - } - - /// - /// Is even? - /// - /// Number to check - /// true if even, otherwise false - static bool IsEven(int number) - { - return number%2 == 0; - } - - /// - /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the - /// solution matrix and X is the unknown matrix. - /// - /// The coefficient matrix, A. - /// The solution matrix, B. - /// The result matrix, X. - public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) - { - var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); - Solve(matrix, input, result, iterator, preconditioner); - return result; - } - /// /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the /// solution matrix and X is the unknown matrix. @@ -408,5 +322,33 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Solvers } } } + + /// + /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the + /// solution vector and x is the unknown vector. + /// + /// The coefficient matrix, A. + /// The solution vector, b. + /// The result vector, x. + public Vector Solve(Matrix matrix, Vector vector, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = new DenseVector(matrix.RowCount); + Solve(matrix, vector, result, iterator, preconditioner); + return result; + } + + /// + /// Solves the matrix equation AX = B, where A is the coefficient matrix, B is the + /// solution matrix and X is the unknown matrix. + /// + /// The coefficient matrix, A. + /// The solution matrix, B. + /// The result matrix, X. + public Matrix Solve(Matrix matrix, Matrix input, Iterator iterator = null, IPreconditioner preconditioner = null) + { + var result = matrix.CreateMatrix(input.RowCount, input.ColumnCount); + Solve(matrix, input, result, iterator, preconditioner); + return result; + } } } diff --git a/src/Numerics/LinearAlgebra/Solvers/IIterativeSolver.cs b/src/Numerics/LinearAlgebra/Solvers/IIterativeSolver.cs index f9b5640f..066b1fd4 100644 --- a/src/Numerics/LinearAlgebra/Solvers/IIterativeSolver.cs +++ b/src/Numerics/LinearAlgebra/Solvers/IIterativeSolver.cs @@ -38,14 +38,6 @@ namespace MathNet.Numerics.LinearAlgebra.Solvers /// public interface IIterativeSolver where T : struct, IEquatable, IFormattable { - /// - /// Stops the solve process. - /// - /// - /// Note that it may take an indetermined amount of time for the solver to actually stop the process. - /// - void StopSolve(); - /// /// Solves the matrix equation Ax = b, where A is the coefficient matrix, b is the /// solution vector and x is the unknown vector.