// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // // Copyright (c) 2009-2014 Math.NET // // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without // restriction, including without limitation the rights to use, // copy, modify, merge, publish, distribute, sublicense, and/or sell // copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following // conditions: // // The above copyright notice and this permission notice shall be // included in all copies or substantial portions of the Software. // // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT // HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, // WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR // OTHER DEALINGS IN THE SOFTWARE. // using System; using System.Diagnostics; namespace IStation.Numerics.LinearAlgebra.Solvers { /// /// Defines an that monitors residuals as stop criterion. /// public sealed class ResidualStopCriterion : IIterationStopCriterion where T : struct, IEquatable, IFormattable { /// /// The maximum value for the residual below which the calculation is considered converged. /// double _maximum; /// /// The minimum number of iterations for which the residual has to be below the maximum before /// the calculation is considered converged. /// int _minimumIterationsBelowMaximum; /// /// The status of the calculation /// IterationStatus _status = IterationStatus.Continue; /// /// The number of iterations since the residuals got below the maximum. /// int _iterationCount; /// /// The iteration number of the last iteration. /// int _lastIteration = -1; /// /// Initializes a new instance of the class with the specified /// maximum residual and minimum number of iterations. /// /// /// The maximum value for the residual below which the calculation is considered converged. /// /// /// The minimum number of iterations for which the residual has to be below the maximum before /// the calculation is considered converged. /// public ResidualStopCriterion(double maximum, int minimumIterationsBelowMaximum = 0) { if (maximum < 0) { throw new ArgumentOutOfRangeException(nameof(maximum)); } if (minimumIterationsBelowMaximum < 0) { throw new ArgumentOutOfRangeException(nameof(minimumIterationsBelowMaximum)); } _maximum = maximum; _minimumIterationsBelowMaximum = minimumIterationsBelowMaximum; } /// /// Gets or sets the maximum value for the residual below which the calculation is considered /// converged. /// /// Thrown if the Maximum is set to a negative value. public double Maximum { [DebuggerStepThrough] get => _maximum; [DebuggerStepThrough] set { if (value < 0) { throw new ArgumentOutOfRangeException(nameof(value)); } _maximum = value; } } /// /// Gets or sets the minimum number of iterations for which the residual has to be /// below the maximum before the calculation is considered converged. /// /// Thrown if the BelowMaximumFor is set to a value less than 1. public int MinimumIterationsBelowMaximum { [DebuggerStepThrough] get => _minimumIterationsBelowMaximum; [DebuggerStepThrough] set { if (value < 0) { throw new ArgumentOutOfRangeException(nameof(value)); } _minimumIterationsBelowMaximum = value; } } /// /// Determines the status of the iterative calculation based on the stop criteria stored /// by the current . Result is set into Status field. /// /// The number of iterations that have passed so far. /// The vector containing the current solution values. /// The right hand side vector. /// The vector containing the current residual vectors. /// /// The individual stop criteria may internally track the progress of the calculation based /// on the invocation of this method. Therefore this method should only be called if the /// calculation has moved forwards at least one step. /// public IterationStatus DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector) { if (iterationNumber < 0) { throw new ArgumentOutOfRangeException(nameof(iterationNumber)); } if (solutionVector.Count != sourceVector.Count) { throw new ArgumentException("All vectors must have the same dimensionality.", nameof(sourceVector)); } if (solutionVector.Count != residualVector.Count) { throw new ArgumentException("All vectors must have the same dimensionality.", nameof(residualVector)); } // Store the infinity norms of both the solution and residual vectors // These values will be used to calculate the relative drop in residuals // later on. var residualNorm = residualVector.InfinityNorm(); // This is criterion 1 from Templates for the solution of linear systems. // The problem with this criterion is that it's not limiting enough. For now // we won't use it. Later on we might get back to it. // return mMaximumResidual * (System.Math.Abs(mMatrixNorm) * System.Math.Abs(solutionNorm) + System.Math.Abs(mVectorNorm)); // For now use criterion 2 from Templates for the solution of linear systems. See page 60. // Check the residuals by calculating: // ||r_i|| <= stop_tol * ||b|| var stopCriterion = _maximum*sourceVector.InfinityNorm(); // First check that we have real numbers not NaN's. // NaN's can occur when the iterative process diverges so we // stop if that is the case. if (double.IsNaN(stopCriterion) || double.IsNaN(residualNorm)) { _iterationCount = 0; _status = IterationStatus.Diverged; return _status; } // ||r_i|| <= stop_tol * ||b|| // Stop the calculation if it's clearly smaller than the tolerance if (residualNorm < stopCriterion) { if (_lastIteration <= iterationNumber) { _iterationCount = iterationNumber - _lastIteration; _status = _iterationCount >= _minimumIterationsBelowMaximum ? IterationStatus.Converged : IterationStatus.Continue; } } else { _iterationCount = 0; _status = IterationStatus.Continue; } _lastIteration = iterationNumber; return _status; } /// /// Gets the current calculation status. /// public IterationStatus Status { [DebuggerStepThrough] get => _status; } /// /// Resets the to the pre-calculation state. /// public void Reset() { _status = IterationStatus.Continue; _iterationCount = 0; _lastIteration = -1; } /// /// Clones the current and its settings. /// /// A new instance of the class. public IIterationStopCriterion Clone() { return new ResidualStopCriterion(_maximum, _minimumIterationsBelowMaximum); } } }